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ANALYSIS AND TESTING OF HIGH ENTRAINMENT 
SINGLE-NOZZLE JET PUMPS WITH VARIABLE-AREA MIXING TUBES 

By Kenneth E. Hickman, Philip G. Hill, and Gerald B. Gilbert 


SUMMARY 

The use of jet pumps is of increasing interest for boundary layer con- 
tiol or control force augmentation in V/STOL aircraft. In typical applications, a 
small mass flow of primary air at pressures up to 400 psia can be used to entrain 
a much larger mass flow of secondary air at ambient conditions. The primary 
nozzle flow is supersonic while the secondary flow is subsonic. The jet pump 
system design objectives may be maximum entrainment, maximum thrust aug- 
mentation, or some combination of the two. Little information is available in the 
literature to guide the designer of jet pumps for such applications. 

In this investigation, an analytical model was developed to predict the 
performance characteristics of axisymmetric single-nozzle jet pumps with variable 
area mixing tubes. The primary flow may be subsonic or supersonic. In the region 
upstream of the section where the central jet reaches the wall, the analysis is based 
upon the hypothesis that the mixing phenomenon is fundamentally similar to the 
mixing of a free turbulent jet with the surrounding fluid. The eddy viscosity values 
used in the analysis are adjusted to allow for the effect of the duct walls on the 
mixing process. Integral techniques are employed in a computer program to solve 
the continuity, momentum, and energy equations to determine the variation of flow 
properties along the mixing tube. Wall boundary layer effects are included in the 
analysis. 


Downstream of the section where the jet reaches the wall, the velocity 
profile is assumed to approach asymptotically the profile for fully developed turbulent 
flow in a pipe. Viscous forces are present throughout the flow so no distinct boundary 
layer analysis is employed. The eddy viscosity is assumed to approach the fully- 
developed flow value asymptotically. Wall friction forces are calculated from the 
fully-developed pipe flow friction coefficient. Integral techniques are employed as 
before to determine the variation of flow properties along the mixing tube. 



An experimental program was conducted to measure mixing tube wall 
static pressure variations, velocity profiles, and temperature profiles in a variable 
area mixing tube with a supersonic (M = 2.72) primary jet. Static pressure varia- 
tions were measured at four different secondary flow rates. These test results 
were used to evaluate the analytical model. 

The analytical model yields good predictions of wall static pressure 
distributions, velocity profiles, and temperature profiles along the mixing tube. 
Therefore, the analysis is believed to be ready for use to relate jet pump perfor- 
mance characteristics to mixing tube design. 
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Section 1 


INTRODUCTION 


1. 1 Background 

A number of STOL aircraft boundary layer control systems now under 
consideration employ jet pumps to entrain large flows of secondary air and direct 
them over deflected flaps to achieve lift augmentation. Some proposed VTOL air- 
craft systems also employ jet pumps for direct lift or control force augmentation. 
The primary, high-pressure flow for the jet pumps can be provided either by bleed 
from main engine compressors or by an auxiliary power unit. 

The use of jet pumps as primary components of V/STOL aircraft 
systems makes necessary the development of new design techniques for these 
devices. In aircraft applications, it is essential to be able to minimize the size 
of jet pumps for particular primary and secondary flow conditions. Jet pumps 
for boundary layer control systems generally must have high entrainment ratios— 
the secondary flow often must be over 10 times larger than the primary flow — but 
pressure rises of only a few psf are needed. The primary flow may be highly 
supersonic. Thus, design procedures which have been developed in the past for 
the more conventional low-entrainment, high-pressure rise industrial jet pumps 
are not suitable for V/STOL aircraft jet pump design. 

1.2 Previous Work 

In an earlier program, Dynatech R/D Company carried out an analy- 
tical and experimental investigation of high- entrainment ratio air-to-air jet pumps 
for the Ames Research Center of NASA (reference 1). This investigation was 
limited to jet pumps with constant-diameter mixing tubes. An analytical pro- 
cedure and a computer program were developed to predict the performance of 
such a jet pump over a range of operating conditions. The accuracy of the analysis 
was confirmed by comparing predicted performance to test results for a number 



of multiple-nozzle jet pump configurations at different primary flow pressure and 
temperature levels. Procedures were demonstrated for matching a jet pump to 
its duct system for maximum entrainment or thrust augmentation. 

The selection of the constant- diameter mixing tube configuration allowed 
considerable simplification of the analysis, design, and construction of the jetpump. 
However, it is unlikely that constant-area jet pumps give the best performance for 
all applications. Almost no information is available to indicate the extent of per- 
formance improvements which can be achieved with other mixing tube configurations. 

A method for predicting the flow behavior in jet pumps with arbitrary 
mixing tube shapes and incompressible flows was reported by P.G. Hill (reference 
2). The method is based upon the hypothesis that the mixing phenomenon in a jet 
pump has a fundamental similarity to the mixing of a free turbulent jet with the 
surrounding fluid. Therefore, as in a free jet, the turbulent Reynolds number- 

jet velocity x duct radius 
^ e T 3 eddy viscosity 

—will remain constant with distance as mixing occurs. This is a rather gross 
simplifying assumption but the resulting flow predictions are good. Static pres- 
sure variations and velocity profiles computed on this basis agreed well with test 
data for Helmbold's converging-diverging mixing tube. Once the static pressure 
distribution is known, the jet pump performance can be predicted without further 

difficulty. 

The analytical methods of reference 2 are limited in application to in- 
compressible flow in axisymmetric jet pumps having a single primary jet. These 
analytical methods must be modified to include compressible flow effects if the 
methods are to be useful for the designer of V/STOL aircraft jet pump systems. 
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1.3 


Objectives of This Investigation 


The specific objectives of this investigation are as follows: 

• to develop an analytical procedure for predicting 
the performance of high-entrainment-ratio com- 
pressible flow jet pumps with arbitrary mixing 
tube geometry. 

• to obtain test results with jet pumps having 
variable-area mixing tubes so that the analytical 
methods can be checked. 

The analytical procedure is formulated to allow prediction of the per- 
formance of a particular jet pump nozzle and mixing tube combination over a range 
of primary and secondary flow conditions. To select the best jet pump design for 
a particular application, the analysis can be used to predict the performance for 
a number of different mixing tube shapes. Comparison of the performance char- 
acteristics will show which geometry is best. The off-design performance of the 
jet pump can be determined by using the same analytical procedures. 
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SYMBOLS 


A 

b 

Cf 

C w 

E 

f 0 

h (n) 

So 

g 2 Cn) 

gg (*?) 

H 

k 

k l 

m 


n s 


P 

P 

Po2 


R 

Ro 


area, ft 2 

u i 

diameter of jet at which U = U G + — ^ — , ft 

wall friction coefficient 

nozzle flow coefficient 

dimensionless eddy viscosity = 775 — 

free jet profile value = f G ( 77 ); equation (1) 

velocity profile at the end of Part 1, equation (7) 

dimensional constant = 32.2, lbm-ft/lbf-sec 2 

auxilary velocity profile, equation (7) 

velocity profile for fully-developed flow in a pipe 

boundary layer shape factor 

specific heat ratio 

suction duct loss coefficient 

entrainment ratio = W 0 /Wj 

number of equal-radial- increment annuli used in integral 
analyses, equation (36) 

static pressure, lbf/ft 2 

stagnation pressure, lbf/ft 2 

secondary flow stagnation pressure after correction 
for suction duct losses, lbf/ft 2 

tube radius, ft 

radius at nozzle exit section, ft 
gas constant x g Q , ft 2 / sec 2 -° R 
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Re m 

Rex 

R e 

®oo’ ®20 

To 

T oj 

Too 

AT„ 


T 

u 

U c 

Ui 


u 


jo 


Ujoo 

Uo 

U r 

V(J) 

W 0 

W x 

W(J,K) 

x 

x core 


Reynolds number based on mean velocity; equation (50) 

turbulent Reynolds number 

momentum thickness Reynolds number 

parameters defined by equations (29) and (30) 

stagnation temperature at any radius in mixing zone, 0 R 

relative stagnation temperature at centerline of jet, 0 R 

stagnation temperature of flow adjacent to the duct, 0 R 

difference between stagnation temperature at any radius 
in the jet and the stagnation temperature of the surrounding 
flow, 0 R 

temperature ratio = T 0 j/T 00 
velocity, ft/ sec 

velocity at centerline of jet, ft/ sec 

velocity at centerline of jet relative to U Q , ft/ sec 

relative velocity at centerline of jet at end of transition 
section, ft/ sec 

relative velocity at centerline of jet at beginning of transition 
section, ft/sec 

velocity of outer stream, ft/ sec 

velocity ratio for transition zone = Uj 0 /Uj 00 

terms in equation (37) 

mass flow rate, secondary flow, lbm/sec 

mass flow rate, primary flow, lbm/sec 

coefficient matrix; equation (37) 

axial position along mixing tube, ft 

length of the transition zone, ft 



y 

radius, ft 

Y(K) 

derivatives in equation (37) 

y 

velocity profile shape parameter 

6 

width of shear layer, ft 

6* 

boundary layer displacement thickness, 

€ 

eddy kinematic viscosity, ft 2 / sec 

tl 

dimensionless radius = y/<5 or y/R 

e 

boundary layer momentum thickness, ft 

A 

velocity ratio U G /Uj 

p 

density, lbm/ft 2 

T 

shear stress, lbf/ft 2 

V 

kinematic viscosity, ft 2 / sec 

Subscripts 

oo 

value at top-hat section 

1 

primary flow 

core 

dimension at end of transition zone 

eff 

effective radius or area of mixing tube 

m 

value at mean area of transition zone 

noz 

primary nozzle exit area 

SD 

suction duct upstream of mixing tube 
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Section 3 


ANALYSIS 


3.1 Purpose 

The purpose of the analysis developed in this section is to predict the 
performance characteristics of compressible flow jet pumps with variable-area 
mixing tubes. The jet pumps may have supersonic or subsonic primary flow 
issuing from a single nozzle located along the axis of an axisymmetric cylindrical 
mixing tube. The secondary flow and the mixed flow downstream must remain 
subsonic. The primary and secondary flows are taken to be the same perfect gas. 

A particular objective of the analysis is to predict the variation in 
static pressure along the length of the mixing tube. Knowledge of this pressure 
variation allows calculation of the thrust augmentation of the jet pump, an es- 
sential parameter for jet pump application studies. 

3.2 General Description of the Analytical Model 

The analysis is based upon the incompressible flow jet pump analytical 
model developed by Dr. P. G. Hill (reference 2). This analytical model, with its 
associated computer program, was modified in the present study to account for 
compressible flow effects. The formulation of the analytical model is described 
in this section. The computer program which is based upon the compressible 
flow model is described in Appendix B of this report. 

The following initial assumptions are made for the analysis: 

1. The primary and secondary flows are the same 

perfect gas. - 

2. No heat is transferred across the wall of the 
jet pump. 



3. The jet pump consists of an axisymmetric, cylindrical, 
variable-area mixing tube with a single primary nozzle 
located along the axis. 

4. The primary and secondary flow conditions and the 
nozzle geometry are assumed to be such that no 
normal shocks or moisture condensation shocks 
occur in the primary flow. 

5. The secondary flow and the combined flows after 
mixing are assumed to remain subsonic throughout 
the mixing tube. 

6. The velocity of the primary jet at the nozzle exit is 
greater than the velocity of the secondary flow. 

7. The static pressure is constant across any section 
perpendicular to the axis of the jet pump. 

Dr. Hill's analysis identifies three distinct flow regimes in a jet pump. 

These regimes are shown in figure 1; they may be described as follows: 

Part 1 - A region in which the jet is approximately self- 
preserving and is immersed in a potential outer stream which 
may be accelerating or decelerating, depending on the shape 
of the duct and the rate of entrainment of mass into the jet. 

Tfonirnnlatinn 7, one - A possible region in which recirculation occurs, 
following a deceleration of the outer stream. At the beginning 
of this zone the "edge" of the jet has not yet diffused to the wall 
and the secondary fluid recirculates through the jet. The pres- 
sure gradient is generally observed to be negligible in this zone. 
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Part 2 - The region downstream of the point (fairly distinct 
in many cases) at which the jet attaches to the wall. An adverse 
pressure gradient is generally established but the relatively high 
shearing forces near the wall tend to accelerate the fluid against 
the pressure gradient. If there is a zone of recirculation, it is 
terminated in a short axial distance by these high shearing forces. 

In addition to these three regions, there is a relatively short transition 
zone between the nozzle exit and the section at which a subsequently self-preserving 
velocity profile is attained. 


In Part (1) the jet velocity profile can be approximated well 
by 

~ u ~ - f o (£) atanyx (1 > 

J 

where U = velocity at radius y 

U Q = outer stream velocity 

IT = jet relative velocity at centerline 
x =s axial position along mixing tube 

y = radius 

5 = width of shear layer (see sketch) 



LLo 

Velocity Profile at Station x 
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The functional relationship f Q (y/5) is determined quantitatively from 
velocity profiles measured in axisymmetric jets discharging into free space. 

f Q (7?) = 1.0004 - 0.01757? - 8.3821T/ 2 + 16.58067? 3 

-12.78777? 4 + 3.6087? 5 (2) 

where 7j = y/6 (Part 1) 

The same relationship holds in the recirculation zone but the axial pressure gradient 
in this region is assumed to be zero. 

The relationship above is used to describe the velocity profile at a parti- 
cular axial station in Part 1 of the mixing tube flow. The continuity, momentum, 
moment-of-momentum, energy, and boundary layer equations are used to determine 
the changes in Uj, U 0 , (5, temperature, and pressure which occur from station to 
station along the mixing tube. To solve these equations, the temperature profile 
must be known so that the density variations across the section can be determined. 
Following Abramovich (reference 3), the stagnation temperature profile is taken to 
be the square-root of the velocity profile. 


where T 

o 


T 


oo 



T o " T oo _ f 1/2 /y\ 

T . o \6 / 

OJ 

stagnation temperature at 
any radius in the mixing zone 

stagnation temperature of sur- 
rounding secondary flow 

relative stagnation temperature at 
center of jet 
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The solution of the moment-of-momentum equation requires shear 
stress values to be known as a function of radius. These values are obtained as 
follows : 

t = e p 3U/ 3y (4) 

where r = shear stress 

e = eddy kinematic viscosity » 

P = density 


The value of the turbulent Reynolds number is assumed to remain constant across 
the flow at any axial station in Part 1. This allows calculation of the eddy viscosity 
from the following equation: 


e = Uj6/Re T 

where Re-p = turbulent Reynolds number 


( 5 ) 


At the beginning of Part 1, the jet mixing process is not significantly affected by 
the presence of the mixing tube walls. Therefore, the value Re-pp = 147, from 
incompressible free jet mixing tests, can be used. Further downstream in Part 1, 
as the jet approaches the walls, the mixing process is altered from a free jet to a 
free wake type of mixing. The change in the mixing process is accounted for by 
using the following equation to determine the eddy viscosity at any station in Part 1: 


[ x + I a-®' 1 - 1 **] w 

IT 

where > - 2 — 

u i 

Boundary layer growth must be taken into account in order to predict 
wall static pressure variations with accuracy. Boundary layer displacement 
thickness variations are obtained in the analysis by using the methods of Moses 
(reference 4). The equations used are described in Section 3.4 in this report. 



In Part 2, the jet has reached the wall. The free jet mixing velocity 
profile is no longer appropriate. Instead, the velocity profile is assumed to follow 

the relationship: 

U/U c = f2 (V) + 7 g2 ( r ?) ^ 

where U c = jet velocity at centerline 

f 2 ( i ) ) = velocity profile at the end of Part 1 

V = y/R 

R = mixing tube radius at the axial position considered 
y = y (x) adjustable shape parameter 
g 2 ( 77 ) = auxilary velocity profile 

At the beginning of Part 2, y is set equal to zero and the velocity profile matches 
the velocity profile at the end of Part 1. The auxilary profile g 2 (V) is chosen so 
that, as y approaches 1.0, the U/U c velocity profile approaches the profile for 
fully- developed flow in a pipe. 

g2 (V) = gg(V) ~ £2 (*») 

where ggC 1 ?) = velocity profile for fully-developed flow in a pipe 

No boundary layer calculations are made in Part 2. Viscous forces are 
present throughout the flow so no distinct boundary layer exists. Wall friction forces 
are calculated from turbulent pipe flow correlations. 


The continuity, momentum, moment-of-momentum, and energy equations 
are used to determine the changes in U c , y, temperature, and pressure which occur 
with distance along the mixing tube in Part 2. The solution of the moment-of- 
momentum equation requires determination of the eddy viscosity as a function of 
radius and axial position. Because the flow in Part 2 becomes asymptotic to fully- 
developed pipe flow, the eddy viscosity must be asymptotic to the fully-developed 
flow value. 
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( 9 ) 


r /t wa ii = y/R = V as y (x) approaches 1.0 


E 2 f - 


L C fd * 

2 9 

~drf ZZ W 


( 10 ) 


where 


E 2 f = C2f/U C R = dimensionless eddy viscosity distribution 
e 2 j = eddy kinematic viscosity for fully-developed pipe flow 
Cfd = T wa H - wall friction coefficient 


An arbitrary function is used to make the eddy viscosity distribution in Part 2 con- 
tinuous with that at the end of Part 1. 

% - Ei (l - y 2 ) + y 2 E 2 f (11) 

where Ef = dimensionless eddy viscosity at the end of Part 1 
calculated from equation (6) 

The paragraphs above have described the basic approaches used for the 
analysis of flow behavior in the variable-area compressible flow jet pump. The 
fundamental assumptions for the analysis have been identified. The sections which 
follow present the sets of equations which must be solved in each of the three regions 
of the flow; the transition zone, the region upstream of jet attachment to the wall 
(Part 1), and the region downstream of the point of attachment (Part 2). 

3. 3 Transition Zone Analysis 

The transition zone begins at the primary nozzle exit plane and has a 
length of approximately 20 jet nozzle diameters. At the nozzle exit plane, the 
static pressure in the supersonic primary flow may be different from the static 
pressure in the surrounding secondary flow. We assume that before mixing of 
the two flows begins, the primary jet expands or contracts isentropically until 
its static pressure matches that of the secondary flow. At the station where this 
accommodation is complete, the velocity profile is assumed to resemble a "top- 
hat" as shown in figure 2. Then mixing of the primary and secondary flows begins. 
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The transition zone continues downstream to the section where the potential core 
in the jet ends. At this point, the f 0 (y/S) profile has been attained and the stagnation 
pressure at the center of the jet begins to fall because of mixing with the secondary 
flow. 


The flow conditions at the end of the transition zone are determined 
by solving three simultaneous non-linear algebraic equations which are developed 
from the continuity and momentum equations written for the transition zone as a 
control volume, and from the condition that the stagnation pressure remains con- 
stant along the centerline of the primary jet. 

The length of the transition zone is measured from the primary nozzle 
exit section to the point where the f 0 (y/ 5) profile is attained. This length is 
designated as x core and must be specified as input data for the analysis. For in- 
compressible flow, equation (12) may be used (reference 3). 

x C ore = 4.08 5 0 (1 + ~^~ ) (12) 

u joo 

where <5 0 = radius of primary jet at top-hat section 

U 00 = secondary flow velocity at top-hat section 

Ujoo = primary jet relative velocity 

For compressible flow with a supersonic primary jet, the value of x core will 
depend on whether the jet is under- or over -expanded as it leaves the nozzle. A 
suitable replacement for equation (12) is not known to be available, so x core was 
arbitrarily chosen to be equal to the mixing tube inlet diameter. This length is 
equivalent to about 18 primary nozzle diameters. 

The transition from the top-hat profile to the fo (y/ 5) profile is assumed 
to occur in a control volume of essentially constant area. The effective mixing tube 
radius at x CO re is calculated by taking the boundary layer thickness into account. 


16 



(13) 


Reff Rcore “ ®^o 


where 0=9 O + 0.001x core (14) 

R e ff = effective radius of mixing tube at x core 
R CO re = ra dius of mixing tube at x core 
Q = boundary layer momentum thickness at x core 
0 Q = inlet boundary layer momentum thickness 
H 0 = inlet boundary layer shape factor = 1.4 assumed 


The flow area available for the secondary flow at the top-hat section is given by 
equation (15). 


A e ff - 7r R 2 e ff - A noz A n oz - A primary flow (15) 

where Aeff = secondary flow area at top-hat section 
Anoz = area of primary nozzle exit section 

The velocity of the secondary flow at the top-hat section is calculated from equation 
(16). 


Uoo = 




Pn A 


eff 


(16) 


where 



secondary flow velocity at top-hat section 
mass flow rate of secondary flow 
density of secondary flow 


The value of p Q in equation (16) is the density corresponding to the local static pres- 
sure and temperature. It is computed by an iterative process using the known values 
of inlet stagnation pressure and temperature and the appropriate perfect gas relation- 
ships. The same calculation yields the value of the local static pressure. 
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follows : 


The primary flow conditions at the top-hat section are calculated as 


T 


1 


ol 



( 17 ) 


where 



static temperature in primary flow at top-hat section 
specified primary flow stagnation temperature 
specified primary flow stagnation pressure 
static pressure from secondary flow calculations 


where 


where U. 

joo 


u i - V 2 TTT R g< T ol - T i> 

primary flow velocity at top-hat section 


primary jet relative velocity 


( 18 ) 


( 19 ) 


The flow conditions at the end of the transition zone are computed by 
using the continuity and momentum relationships and the assumption that the stag- 
nation pressure is unchanged at the center of the jet. The stagnation pressure of the 
secondary flow outside the mixing region is assumed to remain constant during 
transition. The stagnation temperature of the secondary flow outside the mixing 
region, and the stagnation temperature at the center of the primary jet, are assumed 
to remain constant. 


The values of U 00 . U joo , Pi,W 0 , and Wi are known to begin the analysis 
which determines the velocity profile at the end of the transition zone. The continuity, 
momentum, and constant centerline pressure equations at the end of the zone may 
be written as follows on the next page. 
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R eff 

Continuity: 2n j pUydy = W q + 


( 20 ) 


R 


eff 


Momentum 


w,u, w u 

11. o oo 


: 27t f J*L- y dy + (p -P JA - (p,-P ) A + — 2 + — r 

J g oo m '*1 oo' m g„ g, 


( 21 ) 


Constant Stagnation Pressure Along Centerline: P ^ = constant 


( 22 ) 


where A m = 7rR 2 e ff 

The velocity profile at the end of the transition zone is given by equations (1) and 
(2). The temperature profile is given by equation (3). 


To permit equations (20), (21), and (22) to be solved simultaneously 
using standard computer subroutines, these equations were rewritten in terms of 
the following dependent parameters : 


U r = 


U J^ 


u 


joo 


(23) 


X = 


U, 


u 


jo 


(24) 


6 

R eff 



(25) 
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The continuity, momentum, and constant centerline stagnation pressure 
equations in final form are as follows: 

Continuity: 


Let 


'mass 


\V + w, 

o 1 

2 

7 tR , , . p U . 
eff oo 300 


(26) 


It can be shown that 

Cmass = [ u r <|) 2 ( Z 1 - S 20 > + u r S 2 ol 

P OO L ^ J 

k 


where 


= (1 - Soo u?x 2 ) 

00 1 


k-1 


®oo 


k-1 


U' 


joo 


kR T 
g 00 


S 20 1 - S‘W 

00 r 


Z 


1 


1 (X + f Q ) 27 jd 7 ) 

J 1 + Tf 0 V2 - S 00 U^(A + tf 


(27) 

(28) 

( 29 ) 

(30) 

(31) 


IT = 



T 


00 


Momentum: 


Let 


«o<P- p „o» Am +W l u l + w 0 U o 


mom 


7 T R r , p U. 
eff K oo 300 


(32) 


It can be shown that 


mom 


= ( P 


1) 


OO 


00 


R T 
g OQ 

u 2 

300 


(33) 


p ?<l- 


) 2 (Z 2 - as 20 ) + u 2 xs. 
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20 



where 


Z 2 = 


1 2 

(X + f ) 2 V drj 


-J 


+ Tf 1/2 -S U 2 (A+f ) 2 

o oo r ' o' 


(34) 


Constant Stagnation Pressure: 


2 , 2 
1 - S U A 
oo r 


ol 


k -1 


1 - So ° T ° 0 U 2 (1 + A ) 2 


P -AP cr . 
oo SD 


ol 


(35) 


where AP„,„ = suction duct losses (see section 4.2.3) 

oD 


The Z integrals are evaluated by using the following summation: 


n s N. 

z ‘ = VE 't' k 2 ”‘ 


i = 1 


(36) 


where n = number of summation strips, each of the same A 17 
s 

N., D. are defined as follows: 
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Equations (27), (33), and (35) are solved simultaneously to yield values of U , A, 
and 6 /R at the end of the transition zone. The value of the static pressure at the 
end of the zone is then determined from equation (28). 



3.4 Flow Analysis Upstream of Jet Attachment 


In Part 1, the zone between the end of the transition zone and the 
section where the jet reaches the wall, seven variables are determined by in- 
tegral techniques. These dependent variables are U., X = U /U., <5, the static 

J J 

pressure p, the relative stagnation temperature at the jet centerline T^, the 
boundary layer momentum thickness 0, and the boundary layer shape factor H. 

The values of these variables are obtained by solving seven simul- 
taneous equations of the following general form. 


where W(J,K) 
Y(K> 

V(J) 


7 

£ W ( J, K) x Y(K) =V(J) 

K = 1 

a coefficient matrix 

the derivatives of the dependent variables with 
respect to^/R Q 

a set of terms not containing any of the dependent 
variables, evaluated at the x/R q station of interest 
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The Y(K) values are listed below: 



The J simultaneous equations used to evaluate these derivatives 

j = 1 : P QO = stagnation pressure in the flow outside the jet = 

2 : momentum equation for the complete flow 

3 : continuity equation 

4 : energy equation 

5 : moment-of-momentum equation 

6 : boundary layer momentum equation 

7 : boundary layer moment-of-momentum equation 


(38) 


are as follows: 
constant v 


> (39) 


These equations and the W(J,K) coefficients are given in detail in Appendix A. 



The velocity profile for the jet in Part 1 is given by equations (1) and (2), 
with the distribution function f taken from free jet data (reference 5). The jet 
temperature profile is given by equation (3). In the jet (o < r < 6) the shear is 
obtained from equation (4) with the eddy viscosity given by equation (6). 


Outside the jet (6 < y < R), wall shear forces are assumed to be 
negligible in the momentum equation for the complete flow. The boundary layer 
momentum thickness 9 and the shape factor H are calculated from the following 
equations : 



(2 + H) 



dU 

o 

dx 



(40) 


dH 

dx 


-H(H + 1 ) (H - 1) 1 

2 U 


dU 

o 

dx 


h 2 -i 

0 


TTP 

f 0.06 (H-l) 

2 (H + SJRq 0 - 1 


( 41 ) 


The friction coefficient in these equations is taken from the Ludwieg- Tillman 
skin friction equation (reference 6) : 


C 


f 


0. 246 R °- 268 10 


-. 678H 


(42) 


where R fi - Reynolds number based on momentum thickness 


These equations are based on the assumption that the outer (potential) flow at 
velocity U Q is incompressible, and use has been made of the relation between 
the boundary layer and jet mixing parameters given in equation (43) : 



J_ _dA_ J_ 

A dx U. 


dU. 

1 

dx 


(43) 
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The equations above allow the boundary layer development to be 
calculated simultaneously with the main flow mixing. Thus, the boundary layer 
displacement thickness is taken into account when the momentum, continuity, 
and energy equations are integrated across the mixing tube cross section. 

The seven equations (39) are solved simultaneously to yield the values 

of the derivatives (38). Then the derivatives are integrated using Runge-Kutta- 

Merson techniques. This integration yields the desired values of U., U Q , 6, p, 

T ., 9, and H at selected values of x/R along the mixing tube, 
oj o 

If a region of recirculation is present, the value of U Q becomes nega- 
tive. The development of an analysis for the flow behavior in a recirculation 
zone was not included in this investigation. 


3. 5 Flow Analysis Downstream of Jet Attachment (Part 2) 


After the jet reaches the wall, the jet velocity profile is assumed to 
follow the relationship given in equation (7). At the beginning of Part 2, the value 
of the shape parameter y(x) is set equal to zero and the velocity profile is given 
by f 2 ( 77 ). The functional relationship f2 is defined so as to be identical to the final 
velocity distribution in Part 1. 


*0 cn) + xi h£ 01) 

*2(1)= — m 


(44) 


In this equation, f^j? ( 77 ) is the boundary layer profile at the end of Part 1, estimated 
by using a power law for the boundary layer: 



The exponent n must satisfy the values of momentum and displacement 
thickness calculated for the boundary layer at the end of Part 1. The resulting 
equation for follows. 
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1 - rj 


(46) 


where 


R (Hj - 1) 
9, Hj (Hj + 1) 



% = boundary layer momentum thickness at the end of Part 1 
H-^ = boundary layer shape factor at the end of Part 1 


The velocity profile in Part 2 includes the auxilary profile g 2 (??). This 
profile is defined by equation (8) so that, as y approaches 1. 0, the Part 2 velocity 
profile asympototically approaches the profile for fully-developed turbulent flow in 
a pipe. 


In Part 2, the values of five variables are determined by integral 
techniques. These dependent variables are the velocity profile values U c and y , 
the static pressure p, the temperature ratio T, and the temperature of the flow at 
the outer radius of the mixing tube, T 00 . The values of these variables are obtained 
by solving six simultaneous equations of the general form given by equation (37). 

The dependent Y(K) variables are listed below: 


Y(l) = 


a / Uc 

9 <urr> 


jo 
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9 (- 


Ro 
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Y(5) = 


9 (■ 


Toj 
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(47) 


Y(3) = 


9 y 


9(- 


Ro 


Y(6) 


Q / Too 4 
9 (t* — r) 


9(- 


OOl 


Ro 


where Pooi an( 3 Tooi the constant values of Poo a ftd Too in Part 1. The variable 
Y(2) above remains zero throughout the Part 2 analysis; this variable is a redundant 
parameter which remains from an earlier version of the computer program. 
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The equations used to evaluate these derivatives are as follows: 


J = 1 
J= 2 
J = 3 
J = 4 
J = 5 
J = 6 


continuity equation 
energy equation 

momentum equation for the complete flow 

(48) 

moment-of-momentum integral equation 
centerline velocity-temperature relationship 
wall velocity = 0 


These equations and the W(J, K) coefficients are given in detail in Appendix A. 


The form of the stagnation temperature profile must be known in order 
to solve the first four equations (48). For simplicity, the temperature profile was 
assumed to be the same as in a free jet. 


■ T V T° - ° - ■ = V f Q (?t 7 in Part 2 

A oj 


(49) 


This approximation is justified by the test results in Section 4.3 of this report. 


Wall shear stresses are included in the momentum equation for the 
complete flow. The wall friction coefficient used in the analysis is based upon pipe 
flow correlations which yield equation (50). 


Cfdf 


where Re m = 
U 

u c = 

V 


T wall 


U .2 „ -.20 


= 0. 048 ( ) Re 

Ur 


m 


U D 


Reynold’s number based on mean velocity 


mass-average mean velocity 
centerline velocity 
kinematic viscosity 


(50) 
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This wall friction coefficient is only an approximate ~n tothe actual value because 
the velocity profile near the wall in Part 2 oflhe jet pump mixing tube is generally 
not identical to the fully-developed pipe flow velocity profile. Comparison of 
analytical predictions to measured wall static pressure values indicated that 
equation (50) gave values of C fd which were too high. Therefore, the analysis now 

employs an arbitrarily reduced friction factor. 


C 


fd 


1/2 c fd f 


(51) 


The moment-of-momentum integral equation includes a term which 
represents axial shear forces between adjacent stream tubes. These shear forces 
are determined from the eddy viscosity relationship given in equation (11). 


The fifth equation in the set, the centerline velocity-temperature rela- 
tionship, is based upon the test results obtained during this investigation. As shown 
in Section 4. 3 and Figure 16, the following equation may be used to supplement the 
energy equation ip. Part 2. 


j_5l = _J_ 

T dx U dx 

j c 


(52) 


The sixth equation (48) sets Uo, the velocity of the flow along the mixing 
tube surface, equal to zero. This equation was added to eliminate Y(2), the redun- 
dant variable in equation (47), during the solution of the six simultaneous equations 
(48). The solution of these equations yields the values of the derivatives (47). The 
derivatives then are integrated using Runge-Kutta-Merson techniques. This inte- 
gration yields the desired values of Uc, 7, P, Toj, and Too at selected values of 
x/Rq along the mixing tube in the region after the jet reaches the wall. 
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Section 4 


TEST PROGRAM 


The objective of the test program was to provide data which could be 
used to evaluate the analytical model. The test conditions are summarized below: 

Primary Flow 


stagnation pressure: 
stagnation temperature: 
nozzle throat area: 
nozzle geometry: 
mass flow rate: 


348 psia 
807° F 

1.587 x 10~ 4 ft 2 
see figure 3 
6. 76 lbm/min 


Secondary Flow 


inlet stagnation pressure: 
inlet stagnation temperature: 
mixing tube geometry: 
pressure rise: 


laboratory ambient (30. 06" Hg) 

laboratory ambient (92 D F) 

see figure 4 

regulated by discharge 
throttling device 


This section of the report describes the jet pump test arrangement, instrumenta- 
tion and data reduction procedures, and the results which were obtained. 

4. 1 Test Arrangement 

The jet pump test arrangement is shown in figure 5. The primary 
flow was supplied by a 2-stage reciprocating compressor. Electrical heaters 
were used to increase the temperature of the flow up to about 800° F. The 
primary flow was delivered to a single nozzle directed along the axis of the 
mixing tube. 

The momentum of the primary flow entrains a secondary air flow 
from the room into the bellmouth inlet and then into the mixing tube. Here, the 


29 



two streams mix together and the stagnation pressure of the secondary stream is 
increased. The flow from the mixing tube passes through a conical diffuser and 
exhausts to the atmosphere through an adjustable t mottling cone. 

The individual components of the experimental jet pump are described 

below: 

1. Calibrated bellmouth inlet section 

This component consists of a wooden bellmouth, 
metal connecting tube, and fiberglass primary flow 
inlet section. The bellmouth differential pressure 
was calibrated in terms of flow rate by using an 
orifice and blower available in the laboratory. The 
calibrated bellmouth permitted direct measurement 
of secondary mass flow rate for all jet pump tests. 

2. Mixing tube 

The mixing tube geometry was chosen rather 
arbitrarily before the computer program became 
available as a design guide. The basic Helmbold 
mixing tube geometry (reference 7) was selected 
because this geometry has been tested thoroughly 
in the incompressible flow regime. The incom- 
pressible results provide a guide to the flow behavior 
which may be expected in the compressible flow regime. 

The Helmbold mixing tube was scaled down so that 
all dimensions were 0.892 times their original values. 

This scale was selected so that the mixing tube would 
match an existing discharge diffuser and the mixing 
tube throat velocity would remain subsonic for all flow 
rates expected in the test program. The smooth curving 
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profile of the Helmbold tube was approximated with 
cones and cylinders as shown in figure 4 for ease 
of fabrication. 

Discharge diffuser 

A conical diffuser with an area ratio of about 2. 8 
and a total included angle of 7. 1° was added at the end 
of the mixing tube to maximize static pressure recovery 
and allow high entrainment ratios to be achieved. Changes 
in the axial positioning of the throttle cone in the diffuser 
exit produced a variable system resistance so that the jet 
pump could be tested over a range of secondary flow rates. 

Nozzle geometry 

Figure 3 shows the geometry of the converging-diverging 
primary flow nozzle. The area ratio from throat to exit 
section is 3.24, the area ratio corresponding to one- 
dimensional isentropic expansion from 350 psia to 14. 7 
psia. When the jet pump was assembled, the exit plane 
of the nozzle was at x = o where x is defined on the mixing 
tube drawing, figure 4. The mixing tube diameter at the 
nozzle exit plane is 5. 341 in. The nozzle flow coefficient, 
according to the definition below, was measured to be 0.929. 

W 1 

^w W 1 (^3) 

ideal 

where = measured nozzle flow rate at design 
pressure and temperature 

W 1 = isentropic flow rate through nozzle throat 
ideal , , . 

at design pressure and temperature; based 
upon one-dimensional flow assumption 
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4.2 


Instrumentation and Data Reduction Procedures 


4.2.1 Instrumentation 

The instrumentation used to determine the performance of the experi- 
mental jet pump is shown on figure 6 and described in table 1. 

The jet pump inlet bellmouth was calibrated for use as a flowmeter. 
The calibration was accomplished by connecting the bellmouth and the suction 
duct to the inlet of a blower by means of an orifice run and throttling arrange- 
ment. The bellmouth flow equation follows: 

W Q = 229.5 yjp^ Al^ ( 54 ) 

where Ah^ = P b differential pressure, in. HgO gage 

g 

p = inlet density, lbm/ft 

b 

The mixing tube was provided with 21 static pressure taps along its 
length. Four additional static pressure taps were located in the discharge dif- 
fuser. Provision was made for traverse probe measurements at five of the static 
pressure tap sections. The location of all of these taps is given in table 2. The 
exact dimensions were measured at several stations in the mixing tube after its 
construction; these dimensions are also given in table 2. 

The Kiel-temperature probe which was traversed to measure the 
velocity and temperature profiles had a stem diameter of 1/8" . The probe was 
small enough so that probe blockage effects were negligible during the traversing. 

4.2.2 Data Reduction Procedures 

The measured data were used to calculate the following jet pump 
parameters : 
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U vs (X) 

T o - (f ) 
P VS (^) 


- jet pump entrainment ratio 


- velocity profiles 


- temperature profiles 


- mixing tube static pressure variations 


The stagnation pressure and temperature profiles were measured at 
all traverse locations in a plane perpendicular to the axis of the primary flow 
feed pipe (see figure 6). At the station in the mixing tube throat (x/R q = 9. 25), 
a traverse also was made in the plane of the feed pipe to confirm that the flow 
was axisymmetric as desired. 

The wall static pressure and the traverse probe stagnation pressure 
and temperature measurements were used with the appropriate compressible flow 
equations to allow calculation of the velocity profiles at traverse stations 2 through 
6. As a result of a thermocouple failure during the test runs, no temperature data 
were obtained at traverse station 1. Because the temperature profile is required in 
order to calculate the velocity profile, it was necessary to prepare an approximate 
temperature profile for this station. The procedure used is described under Test 
Results in Section 4. 3 of this report. 

4.2.3 Suction Duct Losses 


Results from previous tests of the bellmouth and suction duct assembly 
(reference 1) and static pressure data from the present test program indicate that 
stagnation pressure losses in the suction duct upstream of the mixing tube are of 
the order of 2 in. H^O for the tested secondary flow rates. These losses may be 
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accounted for in the jet pump analysis by using equation (55) to calculate the 
secondary flow stagnation pressure at the primary nozzle exit section in the 


mixing tube: 


o2 


= P - K t 
oo L 


TT 2 

P oo U SD 

2g„ 


or 





2g p A 
6 o^oo 


2 

SD 


(55) 


where 



K, 


oo 


U. 


SD 


SD 


= stagnation pressure at primary nozzle exit plane 

= stagnation pressure at suction duct inlet (laboratory 
ambient) 

= suction duct loss coefficient 

= density corresponding to suction duct inlet 
stagnation state 

= suction duct velocity (assumed uniform) 

= suction duct cross-sectional area 


For the suction duct in the experimental jet pump, the value of is 0. 33. 


4. 3 Test Results 

The jet pump was tested at four values of entrainment ratio, 17. 0, 
19,4, 21. 0, and 23. 6. The corresponding values of primary and secondary 
mass flow rates are given in table 3. The inlet pressures and temperatures 
were constant throughout the test and were as listed at the beginning of this 
Section 4. 
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Wall static pressure values measured along the mixing tube are listed 
for all four entrainment ratios on table 3. These values are plotted in figure 7. 

The operation of the jet pump was reasonably steady (i.e. , wall pressure 
fluctuations were small) when the entrainment ratio was 21. 0. Therefore, this 
condition was selected for the velocity traverse measurements which require long 
periods of steady operation. The velocity profiles for traverse stations 2 through 
6 are shown in figure 8. The associated temperature profiles are shown in figure 9. 

Traverses 4 and 5 were taken at the same station in the constant -area 
throat section of the mixing tube. The axes of the traverse were 90° apart so that 
any departures from axial symmetry in the flow could be detected. The slight 
departures which were observed are due to heating of the secondary flow as it passes 
over the primary nozzle feed pipe upstream of the mixing tube inlet. These departures 
have a negligible effect on jet pump performance and will not interfere with our 
comparison of measured and predicted flow behavior through the jet pump. 

Because of a thermocouple failure during testing, no temperature data 
were obtained at traverse station 1. The stagnation pressure measurements at this 
section cannot be used to determine the velocity profile unless the temperature 
profile is available. An approximate velocity profile for traverse station 1 was 
developed by using the analytically-predicted temperature profile together with the 
measured stagnation pressure values. The resulting velocity profile is given at the 
end of the next section of this report. 

The mass flow rate through the jet pump as determined by the calibrated 
inlet bellmouth was compared to the mass flow rate obtained by integration of the 
velocity profiles for stations 4 and 5. Agreement was within 1 % (149. 8 lbm/min. 
from integration vs. 148. 8 lbm/ min. from the bellmouth). The measured velocity 
profile at station 6 was used for a similar comparison. Integration of this profile 

gave a mass flow rate of 158. 8 lbm/ min. , about 7% greater than the bellmouth 
measurement. 
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The values of and calculated from the measured velocity and 

temperature profiles are plotted in figire to to show how the centerline velocities 
and temperatures vary with distance along the mixing tube. The velocity and temp- 
erature ratios are nearly identical over most of the mixing tube length. 
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Section 5 


COMPARISON OF ANALYTICAL AND EXPERIMENTAL RESULTS 

5. 1 Mixing Tube Wall Static Pressure Variations 

The mixing tube wall static pressure measurements are the most valu- 
able results for evaluating the accuracy of the analytical model for use in jet pump 
design. The prediction of this pressure variation is the primary purpose of this 
investigation because knowledge of this variation permits calculation of the pressure 
force on the mixing tube wall. This force must be known in order to solve the 
momentum equation during jet pump system optimization studies. 

The analytical predictions of mixing tube static pressure variations are 

compared to test results for four entrainment ratios in figures 11, 12, and 13. The 

analyses were carried out with two different values assumed for x , the length 

core & 

of the transition region at the primary nozzle exit, and for two values of the second- 
ary flow rates for each test; the values determined from the test results using the 
bellmouth calibration equation (54), and values 2% lower. A key to the three figures 
follows: 


Figure No. 

x /R 

core o 

Secondary Flow Rates 
from (54), reduced by 2 % 

11 

2.5 

12 

2. 5 

from (54) 

13 

2. 0 

from (54), reduced by 2% 


The mass flow rates given by equation (54) and used to prepare figure 12 cause the 
analytical predictions of static pressure to fall below the measured values in the 
throat section of the mixing tube (^- from 7. 34 to 10. 7). The assumption that the 
secondary flow rates are 2% lower yields better agreement as shown in figures 11 
and 13. The choice of x to be 2. 5 R rather than 2. 0 R causes only a small 
difference in the predicted static pressure levels. The differences are largest in 
the diffuser sections downstream of the mixing tube throat. 
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From these results, it was concluded that further comparisons of 
analytical and experimental results should be based on the assumption that the true 
secondary flow rates are 2 % lower than the flow rates given by equation (54). An 
uncertainty of ± 2% in flow rate is not unreasonable for the bellmouth calibration. 

The 2 % flow correction brings the analytical predictions very close to the experimental 
results except for the static pressures downstream of the mixing tube throat. 

5.2 Velocity and Temperature Profiles 

The variation of predicted centerline velocity with distance along the 
mixing tube is shown in figure 14 for three alternative values of 
and 2. 5. The measured values of centerline velocity at traverse stations 2-6 are 
also plotted in the figure. A value of x core / R Q between 2. 0 and 2. 5 appears to 
make the analytical prediction fit the test data most accurately. 

The variation of predicted centerline stagnation temperature with distance 
along the mixing tube also is shown in figure 14. A value of x core /R Q between 2. 0 
and 2. 5 will make the temperature predictions fit the test data upstream of the 
throat section of the mixing tube. At traverse stations 4, 5, and 6, the measured 
temperature levels fall about 30° F below the predicted centerline stagnation 
temperatures. 

The analytical results (U,. U Q , , f 2> g 2> and v), together with the 
known free jet profile f {'<} ), allow direct comparison of the velocity and tempera- 
ture profiles predicted by the analysis to the velocity and temperature profiles 
measured during the test program. The velocity profiles are compared in figure 15, 
and the temperature profiles are compared in figure 16. The predicted velocity 
profiles agree reasonably well with the measured profiles. The measured and 
predicted temperature profiles agree well for traverse stations 2 and 3, but the 
predicted temperatures near the centerline for stations 4, 5, and 6 are somewhat 
higher than they should be. 

Stagnation pressure measurements only were obtained at traverse 
station 1. These measurements, coupled with the analytical temperature profiles 
predicted for this station, can be used to develop an approximate velocity profile- 
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The procedures used were as follows: 


1. The stagnation pressure data from the traverse probe, together 
with the local static pressure tap reading, were used to determine 
the Mach number and t/ T q ratio at each y/ R position in the 
mixing tube cross-section. This data is given in table 4. 

2. From the analytical solution for x = 2. 5 R , the predicted 

^ coi*g o 

value of p — at the traverse station (■— - = 2.5) was found to be 

0..2118, T$e local value of p is 0. 8$9. These results allow the 

y/ R positions of the traverse probe near the duct centerline to be 

interpreted in terms of the y/ <5 values for the free jet velocity 

profile of equations (1) and (2). Using the analytically predicted 

values U = 3019 ft/ sec, U = 268 ft/ sec, T = 1267° R, and 

T 00 - 552 R, the free jet velocity and temperature profiles can 

be used, through equations (1), (2), and (3), to determine the 

predicted values of velocity and stagnation temperature for each 

y/R position within the jet mixing region. The corresponding 

static temperatures can be determined from the T/ T ratios in 

o 

table 4. The speed of sound is calculated from the static tempera- 
ture. The predicted flow velocities and speed of sound values are 
used to calculate Mach numbers for each y/R position within the 
jet mixing region. These predicted Mach numbers are compared 
to the "measured " Mach numbers in table 4. If the predicted and 
measured numbers agree, the associated velocity and temperature 
profiles afford a good approximation to the true profiles. 

t. The same calculation procedure was followed using the analytical 
solution for x cQre =2.0 R q . The predicted value of at the 
traverse station was 0.287. The other predicted values°employed 
in the analysis were as follows: 

U = 2227 ft/ sec. T = 1041° R 

c oc 

U = 261 ft/ sec. T = 552° R 

° oo 
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The predicted and "measured” Mach number profiles for traverse station 1 are 

compared in figure 17. The predicted profile based upon the assumption that 

v = 2. 0 R is closer to the measured profile than the x = 2. 5 R profile 
A C ore o cuic w 

although the predicted centerline velocity is too high. The predicted velocity and 
temperature profiles for both x CQre assumptions are given in table 4. The 
predicted values were obtained using a secondary flow rate which was 2% less than 
the value given by equation (54). 
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Section 6 
CONCLUSIONS 


An analytical method has been developed to predict the performance 
characteristics of axisymmetric single-nozzle compressible flow jet pumps with 
variable area mixing tubes. The primary flow may be either subsonic or supersonic. 
The analysis is divided into two parts. In part 1, the region between the primary 
nozzle exit and the point where the jet reaches the wall, the analysis is based upon 
the hypothesis that the mixing phenomena in the jet pump is fundamentally similar to 
the mixing of a free turbulent jet with the surrounding fluid. The eddy viscosity is 
adjusted to account for the influence of the duct walls as the jet approaches the walls. 
In part 2, downstream of the point where the jet reaches the wall, the velocity 
profile is allowed to vary from the free jet profile at the end of part 1 to a profile 
which asymptotically approaches the fully-developed turbulent flow profile in a 
pipe. Integral techniques are employed in both part 1 and part 2 to solve the contin- 
uity, momentum, moment-of-momentum, and energy equations to determine the 
variations of flow properties along the mixing tube. 

An experimental program was conducted to measure mixing tube wall 
static pressure variations, velocity profiles, and temperature profiles in a variable 
area mixing tube with a supersonic (M = 2. 72) primary jet. Static pressure varia- 
tions were measured at four different secondary flow rates. These test results were 
used to evaluate the analytical model. 

Analytical predictions of wall static pressure distributions along the 
mixing tube generally agreed well with the test results for all four entrainment 
ratios. The predicted wall static pressure values differed slightly from the 
measured pressures downstream of the constant-area throat section. The velocity 
profiles along the mixing tube were predicted accurately by the analysis. The 
analytical temperature profiles were not as accurate; the predicted centerline 
temperatures downstream of the throat were too high. These discrepancies are 
considered to be minor in view of the comparatively extreme mixing tube geometry 
used for the test case. Thus, the analysis is ready for use to calculate the pressure 
force on the wall of a variable area mixing tube. This permits the momentum 
equation to be solved accurately in jet pump-duct system optimization and design 
studies. 
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The analysis in part 2 of the jet pump makes the assumption that the 
temperature profiles are similar to free jet temperature profiles. A very simple 
and approximate form of the energy equation is employed. A more accurate energy 
equation, perhaps augmented by assumption of a different form for the temperature 
profile, might lead to greater accuracy in the prediction of wall static pressures 
and temperature profiles in this region. 
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APPENDIX A 


Equations for the Flow 


A1 - Part 1 - Upstream of Jet Attachment 

The general form of the flow equations, as described in Section 3. 4, 
is as follows: 


^ W(J,K) * Y(K) = V(J) 
k=l 


The 7 variables are tabulated below, using the convention that the superscript (') 

a 


represents 




K = 1 



Y(K) = X 


1 ( 8 \ ( p V i °i i / e \ 

Ir J Ip 1 It I Ir J 

\o/ \oo/ \ oo / \o/ 


H’ 


The W(J,K) coefficients and V(J) terms are determined in this section. 

Al-1 Equation for J = 1; Constant stagnation pressure in the flow outside the jet. 


dp = - p 


UdU 

S n 


R T-^2- + XU.d(AU.) = 0 

g p i ' r 


Normalizing: 


RT ^ U. 
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= 0 
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Let BP = 


R T 

_£ 00 


k-1 


U 2 2kS o 

jo 


S o " „ k 


UV 


2 k-1 Voo 


Then T o = T oo 


U. 


1 - S 


0 U. 2 

]o 


= Static temperature in the flow 
outside the jet 


The final values follow: 

2 V 

W(l,l)-X 

jo 
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W(l,2) = XI^- 


3° 
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W(l, 6) “ 0 
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Al-2 Equation for J = 2; Momentum equation for the flow 
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where A 
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In the computer analysis , this integration is approximated by a 
summation across the jet: 


Let 




N 


iJ 


D. 


V t 


iJ 


In this equation N. and D. are average values of the numerator and denominator 

i.L 1 1 

across the 1 equal- radius annular segment of the jet. The following additional 

definitions will be used: 
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Then the relations below may be used: 



8Z U 

TT =Z 5J + Z 6J 

Additional parameters which simplify the equations are defined as follows: 
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R _ tube 0 TT 

n r> -r» ” 
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W(2,3) = 2 (Z 12 - XS 2 ) 
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W(2,4) 


W(2, 5) 



W(2,6) =-2^- S 2 X 
W(2,7) = -2-j| S 2 X 

V(2) = - 2 R ^ Ul?g - S 2 X 

Table A1 lists the values of N., D., and their derivatives which are required to 
evaluate the Z parameters in the previous equations. 


Al. 3 Equation for J = 3; Continuity equation 
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The continuity equation is normalized as follows: 
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W(3,7)= -2S 2 ^ 


V(3) = - 2 




Table A1 lists the values of N., D^, and their derivatives which are required to 
evaluate the Z parameters in the equations above. 

Al-4 Equation for J = 4; Energy Equation 
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Then, the normalized energy equation may be written as follows: 
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If this equation is compared to the normalized continuity equation in Section A 1-3, 
it is seen that the right-hand sides are identical except for the substitution of 
Z, . for Z 1<J . This means that all of the W(4,K) coefficients are identical to the 
W(3,K) coefficients except for the substitution of Z.^ for Z ^ in all expressions. 
Table A1 lists the values of N., D., and their derivatives which are required to 
evaluate the Z. 4 parameters. 

Al-5 Equation for J =5; Moment- of- Momentum Integral Equation 

The momentum equation for an annular section of the jet can be derived 


as follows. 
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To derive the moment- of momentum integral equation, this momentum equation is 
multiplied by ydy and integrated across the jet: 


f puy ir ydy + f pvy §7 ydy = f ydy -f ar y2<ly 
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In order to evaluate the radial velocity, v, it is necessary to use the 
continuity relation. Employing a control volume of radius y and length dx, the 
continuity equation may be written as follows: 
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The pressure gradient term may be rewritten as follows: 
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The shear stress term may be evaluated as follows: 



ydy = ry 
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au u i 9f 

t = Pe -gy = P e where e is the eddy kinematic 

viscosity 
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where A = - / -| “FT gpr- ’Id'n = 0. 377 as in incompressible jet mixing 
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fu 6 J * s f nver ® e °f the local turbulent Reynolds number. 
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Assembling all the terms, the final moment- of- momentum integral 
equation is as follows: 
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The values of the coefficients follow: 


U. 


W(5,l) =Q X +R X +R 2 rtf— 

jo 

W(5, 5) 

W(5,2) =Q 2 +R 34 

W(5, 6) 

W(5, 3) =Q 3 +R 1q +R n 

W(5, 7) 

W(5,4) =R X + 

V(5) = 


U/ 
3 3 


Al-6 Equation for J = 6; Boundary Layer Momentum Equation 

The boundary layer momentum equation, equation 
in Section 3. 4 of this report. 



= 0 


= 0 
EA t 

W\) 


(40), is discussed 


Al-7 Equation for J = 7; Boundary Layer Moment- of- Momentum Equation 

The boundary layer moment- of- momentum equation was used to derive 
the shape factor equation (41). This equation is discussed in Section 3.4 of this 
report. 


Al-8 Initial Conditions for the Part 1 Analysis 

The initial conditions for the Part 1 analysis are established from the 
transition zone analysis described in Section 3.3. The initial values were set 
as follows: 
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A2 - Part 2 - Downstream of Jet Attachment 

The general form of the flow equations, as described in Section 3.4, 
is as follows: 


W(J,K) * Y(K) = V(J) 

fel 


The 6 variables employed in Part 2 are tabulated below. The superscript (') represents 

a 


the derivative 


a 




K = 


YW "(v)fi) y ' fcr) T ’ fe) 

where P and T . are the stagnation pressure and stagnation temperature for the 
ooi 001 

wall streamline at the end of Part 1 just as the jet reaches the duct wall. 
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The variable Y(2) remains zero throughout the Part 2 analysis; this 
variable is a redundant parameter which remains from an earlier version of the 
computer program. 

The W(J,K) coefficients and V(J) terms are determined in this section 
of the appendix. 

A2-1 Equation for J = 1 ; Continuity Equation 
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where - T — T ) (Free jet temperature profile) is assumed to 

hold in Part 2 as a simplification of the 
analysis. 


The value of T used in the definitions of BP and S , and throughout the Part 2 
oo o 

analysis, is the stagnation temperature for the wall streamline at the axial position 

selected. T varies with x in Part 2. 
oo 
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Let 
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The continuity equation may be rewritten as follows: 
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In the computer analysis , the integration is approximated by a summation 
across the flow: 
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T) . 
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In this equation, N. and D. are average values of the numerator and denominator 
til 11 

across the i — equal- radius annular segment of the flow. 

The following additional definitions will be used: 
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Employing these definitions in the continuity equation, the following equation is 
obtained after differentiating: 
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Table A2 lists the values of N. , D., and their derivatives which are required to 
evaluate the Z parameters in the previous equations. 

A2. 2 Equation for J = 2; Energy Equation 
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Pu T q • 27 rydy 


assuming constant specific 
heat throughout the flow 


Using the substitutions forp, the velocity profile functions, and q as in Section 
A2. 1, the energy equation may be rewritten as follows: 
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As in Section A2. 1, the free-jet temperature profile is assumed to hold: 


= 1 + 
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With this, the energy equation becomes: 
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Table A2 lists the values of d. , and their derivatives which are required to 
evaluate the Z parameters in the equations above. 

A2. 3 Equation for J =3; Momentum Equation 
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The momentum equation may be rewritten as follows: 
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Collecting terms and dividing each by 





o 


1 

BP : 


W(3 , 1) = 



+ Z 


23 


W(3 , 2) = 0 



W(3 > 5 ) = Z 53 + Z 63 



+ Z 


83 



66 


Table A 2 lists the values of and their derivatives which are required to 

evaluate the Z parameters in the equations above. 



A2. 4 Equation for J - 4; Moment- of- Momentum Integral Equation 


The moment- of-momentum integral equation is taken from Section Al. 5 
of this appendix: 
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In which 



d = i - s a ) <f 2 *rs/ 


Following the analysis In Section Al. 5, the second integral is evaluated as follows: 




68 





69 








70 




in which 
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The pressure gradient term may be written as follows: 
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The shear stress term may be evaluated as follows: 
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where e is the eddy kinematic 
viscosity 



71 



R„ R r e i fn 9 
-(i/ 2 )c fd -R- - -R- J TTr T5 - 


+ 7 


ag, 

9*T 


77d77 


R 


R C.j 

o fd 

R 2 


where A, 


1 


E, 


af 2 9g 2\ 

-^=- + y -z-=- ) rjdr) 
> a^7 _ Or? / 

D 


E f 


U R 
c 


in Part 2 


Assembling all the terms, the final moment-of- momentum integral 
equation is as follows: 
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The values of the coefficients follow: 
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finally. 
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The values of the coefficients follow: 



W(5,l) = 



W(5, 5) 
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W(5, 2) = 0 



W(5, 3) = 0 
W(5, 4) = 0 


V(5) = 0 


A2. 6 Equation for J = 6; Wall Velocity - 0 


The value of U is assumed to be zero throughout Part 2. Therefore, 
o 
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and W(6, 1) = W(6, 3) = W(6,4) = W(6,5) = W(6,6) = V(6) = 0 

W(6, 2) = 1 
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Appendix B 

THE COMPUTER PROGRAM 
B. 1 General Description 

The computer program has 10 sections. The general functions of each 
section are described below. 

MAIN : The program begins and ends in MAIN. Input data concerning the jet 

pump geometry, inlet gas flow properties, free jet velocity profile, and stations 
along the mixing tube where output values are desired, are all read in by MAIN 
and by two subroutines called by MAIN — DIEFEQ and SUB. MAIN converts the 
units of the input parameters into other units which are more convenient for 
subsequent analysis. 

After conversion of units, MAIN computes the primary and secondary 
flow conditions at the top-hat section as described in Section 3. 3 of this report. 
Then MAIN sets up the initial trial values for the velocity profile after transition 
and calls VB04A to perform the iterations required to obtain an accurate solution 
for the profile. 

After the transition zone analysis has been completed, MAIN sets up 
the initial conditions for the flow analysis upstream of the point of jet attach- 
ment to the wall. It also defines the stations along the mixing tube for which 
data will be printed out. MAIN then calls RUNGE to carry out the solution for 
the remainder of the flow analysis. 

SUB: The first section of SUB, called when J - 3, reads in data on the 

mixing tube geometry — inner diameter vs. length. The diameters are converted 
to radii and all radii and length values are made non-dimensional by dividing by 
R . The second section of SUB, called when J - 1 or 2, finds the duct radius 
and slope at any axial position x specified as an input value to the subroutine. 
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The procedure used is linear interpolation between the nearest upstream and 
downstream radii which were read as input data by the first section of SUB. 

CALXFG: The purpose of CALXFG is to perform the computations required 

to set up the three transition zone equations (27), (33), and (35) for solution 
by VB04A. The three equations and derivatives of each of the three equations 
with respect to the three variables U^, X, and 6^ are computed in CALXFG. 

VB04A, VD02A, and SPNIST : These subroutines are library routines employed 

to solve the three simultaneous non-linear algebraic equations (27), (33), and (35). 

A two-page discussion of these subroutines is included at the end of section B. 3. 

DIFFEQ : The DIFFEQ subroutine is divided into two parts. Part I establishes 

the 7 simultaneous equations (39) which must be solved to determine 
the flow conditions upstream of jet attachment. The equations used 
are outlined in section 3.4 and detailed in appendix A. 1. When the simultaneous 
equations are set up, DIFFEQ calls subroutine SIMQ to solve the equations for 
the values of the 7 derivatives in equation (38)- Then subroutine RUNGE is called 
to integrate the derivatives using Runge-Kutta-Merson techniques. This integra- 
tion yields the values of U., U , 5, p, T ., 0, and H at stations closely spaced 
along the duct. 

Part 2 of DIFFEQ establishes the 6 simultaneous equations (48) which 

must be solved to determine the flow conditions in the mixing tube downstream of 

jet attachment to the wall. The equations used are outlined in section 3. 5 and 

detailed in appendix A. 2. Subroutine SIMQ is called to solve for the 6 derivatives 

in equation (47). Then subroutine RUNGE integrates the derivatives to find the 

values of U , y, p, T ., and T at stations closely spaced along the 
C OJ oo 

duct. 

SIMQ: This is a library subroutine which is called by DIFFEQ to solve 

simultaneous linear equations to find the values of the Y(K) derivatives in equations 
(38) and (47). 
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RUNGE: The RUNGE subroutine performs a Runge-Kutta-Merson integration 

procedure to integrate the derivatives of the Y (K) quantities which are developed 
by DIFFEQ and SIMQ as described above. RUNGE also calls the subroutine PRINT 
to print the desired output values of jet pump flow parameters at each mixing 
tube station (XOUT) which has been specified by input data and equations in MAIN. 

PRINT : This subroutine contains instructions for printing the computer jet 

pump flow parameters at selected stations along the mixing tube downstream of 
the transition zone. 

B. 2 Input Data Format 


The input data to the program must be prepared according to the 
following sequence: 


Card No. 

Parameters 

Format 

1 

NS 

312 

2 

GG(I), 1= 1, NS 

10F5.4 

3 

SDLOSS, ASD 

2F10. 3 

4 

THETA, SHAPE, VISC, RZERO 

4F10.6 

5 

1-CARD MESSAGE identifying solution desired 

8 OH 

6 

DELTAX, XTUBE, TURBNO, NSUB, NGAM, 
XCORE, ANOZ 

3F10.4, 215, 
2F10. 5 

7 

POO, TOO, POl, T01 , AMASS1 , AMASSO, 
AG, RG 

8F7.3 

8 to 8 + I 

X(I), A (I) \ 

Q Q V omit if NSUB = 2 

2F15.4 

8+1+1 


2F15.4 


Cards 1 through 7 are required for each solution desired. The cards from 8 on are 
required to define a new mixing tube geometry for analysis. If the same mixing 
tube geometry is to be used for additional solutions with altered flow conditions, the 
cards from 8 on do not have to be included for these additional solutions. The input 
parameter NSUB tells the computer whether the cards from 8 on are included with 
a data set, i.e. , whether the same mixing tube geometry is to be used for additional 
solutions. 
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The input parameters are described below. 


NS 

GG(I) 

SDLOSS 

ASD 

THETA 

SHAPE 

vise 

RZERO 

DELTAX 

XTUBE 

TURBNO 

NSUB 

NGAM 

XCORE 

ANOZ 

POO 

TOO 

POl 

TOl 


number (= 10) of equal-radial-increment strips used to approximate the 
jet mass flow, momentum, and energy integrals across the jet 

average values of -jj- (t] ) for I = NS equal-radial-increment strips for 

the turbulent pipe flow velocity profile 

suction duct loss coefficient; K L in equation (55) 


suction duct area, ft ; 


Agcin equation (55) 


boundary layer momentum thickness at x = 0, ft 
boundary layer shape factor at x = 0 

gas kinematic viscosity for secondary flow at inlet, ft / sec 

mixing tube radius at nozzle exit section; x = 0, ft 

steps of x/ R q at which data printouts are desired in the mixing tube 

mixing tube length, ft 

turbulent Reynolds number value = 147 

control instruction: if 1, a new mixing tube geometry is read in 

if 2, the tube geometry from the previous 
solution will be used. 

control instruction: if 0, incompressible flow solution (not operable) 

if 1 , compressible flow solution 
length of the transition zone divided by R q 

2 

primary nozzle exit flow area, ft 


stagnation pressure upstream of the suction duct losses, psia 


stagnation temperature of the secondary flow, °R 


stagnation pressure of the primary flow, psia 
stagnation temperature of the primary flow, °R 
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AMASS1 mass flow rate of the primary flow, lbm/sec 


AMASSO mass flow rate of the secondary flow, lbm/sec 


AG 


specific heat ratio of the gas 


RG 


gas constant, ft-lbf/lbm-°R 


X(I) 


x stations along the mixing tube at which A(I) values are defined, ft 


A(I) diameter of the mixing tube at the corresponding x , , . , ft 

station 

B. 3 Output Data 

A complete sample of output data from the computer program is given 
in section B. 5 of this appendix. The first section of the output repeats the input 
data and thus requires no comment. The remainder of the data is summarized 
below. 


F (I)’s: values of f Q (77 ) at 77 = 0. 05, 0. 15, 0. 25, 0. 95 

CONDITIONS AT BEGINNING OF THE TRANSITION SECTION 


U 


lists values of U 0O , P o0 , U T 00 , P<)0 (psfa and in.H 2 0), \ , 

joo 

and primary jet momentum = Wj where U^is the velocity achieved 
by isentropic expansion of the primary flow to the static pressure at the 
end of the accomodation process, p . 
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The next portion of the printout monitors the solution by VB04A of 
equations (27), (33), and (35) for the transition zone. Each iteration employing 
CALXFG is recorded. The VARIABLES are the values of U , A , and o/R eff 
determined during the particular iteration reported. The FUNCTIONS are the 
values of the functions: 


mass 


- C 


new 


mass 


old 


mom. 


new 


mom 


old 


const 


- P 


new 


const 


old 


mass 


old 


mom 


old 


const 


old 


where C is defined by equation (27) 
mass 

C is defined by equation (32) 
mom 

k-1 


v oo' 


P .is 
const 


( ) is the value for the current iteration. 

' 'new 

( ) is the value for the previous iteration. 
' old 


If these functions are computed to within ERR times the "old" value of ^ mass > 

C or P , VB04A is judged to have converged satisfactorily. In the present 

mom’ const 

program, ERR is set at 10 , an excessively tight tolerance. As a result, the 

message "VB04A ACCURACY CANNOT BE ACHIEVED" is often printed out. Following 
this message, the values of the VARIABLES and FUNCTIONS for the current iteration 
are printed. These values are used as the first values for subsequent calculations. 


Four lines of print follow the end of the VB04A material. The first 
line restates the values of XX(1) = U y , XX(2) = X, andXX(3)= 6 R gff in numerical 
form. The second line compares the values of EN, a dimensionless jet pump 
parameter developed in reference 2, before and after transition (EN vs. EN2). 


EN 


W + W 
1 o 


V 27TR 2 p 

o oo 


g (p-P )7r R 
6 o' oo' 


eff 


wu, + w u 

11 0 


o 
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The two values should be identical; differences which exist provide a measure of 

the accuracy of the transition analysis. Following the EN values, the values of 

Soand Bp (see appendix A, section Al. 1) are printed. Then the final values of 

U. and U are given, 
jo co & 


The printing continues with a tabulation of values along the mixing tube 
given by Part 1 of the analysis. The parameters listed are as follows: 


X/RZERO; 

AREA; 



PH20; 

UO; 

UCENT; 

UR; 


wall static pressure, in H„0 relative to P 

2 oo 

value of U , secondary flow velocity 
value of U , velocity of flow at the centerline 


value of 


V u 


co 


U 

LAMBDA; value of A = -^p- 

j 

DELTA/R; value of S/R 

o 

TOCENT; value of stagnation temperature at the duct centerline 

TOWALL; value of stagnation temperature in the secondary flow 

outside the mixing region 
THETA/RO; value of 0 /R q 

SHAPE; value of H 


When the Part 1 analysis indicates that the jet reaches the wall, the 
message "DELTA/R = 1 -- DIFFERENTIAL EQUATIONS CHANGE" is printed. 
The local value of U^, called CL, is also printed. Next, two lines are printed as 
follows: 
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F2(I)'s: values of fgfa ) at 77 = 0.05, 0.15, 0.25,- 

G2(I)'s: values of gg(T 7 ) at 77 = 0.05, 0.15, 0.25, 


,0. 95 
-,0.95 


The printing concludes with the tabulated results obtained from the 
Part 2 analysis . The parameters listed are as follows: 


X/RZERO, AREA, PH20, UCENT, UR, TOCENT, TOWALL; same as in Part 1 

TOWALL/TOO; stagnation temperature of wall streamline divided by secondary 
flow inlet stagnation temperature, T 

AUGMENT: the value of the local momentum flux 


00 

R 


f jEL 

’ J So 


27rydy., divided 


by the primary jet momentum, W-^U^, which is printed out earlier. 


GAMMA; local value of y 
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B. 3 Listing 


^^^^*t*i*.****ir*t**i** ************** *************************************** ******* 

C COJET, ANALYSIS OF FLOW BEHAV T fl R IN AXIS YMMETR I C COMPRESSIBLE FLOW 
C EJECTORS WITH VARIABLE ARTA MIXING TUBES 

C* *, * *< ************* <.*******»********* r ***********************♦**+************ * * 

COMMON TURBNO, CF, EN , NP AR T , C , JRUNGE , NG A M 

1, POO, POX , TOOtTOl I AMACH, AG.RG 

2 , AMASS 1, AM AS SO, T 

3 — tSOOir — CMAS, C»CM,- CFNR 

4 ,SO, BP, CM, VISC1 

COMMON ZA, ZB, IJJO, DOR IN ~~ """ ' '~ ' 

6 .MASRAT.RZFRO, AFAC, CFD, HD, F, Al, AAZERO, TF , FP, GP 

7, SDLOSS, ASD 

8, AUG1,U1,UCFNT 

DIMENSION Y ( 10) vTOLUOI ,YMIN(10-r, XOOTt201 1 ,M4<LKT5) 

1 , DY ( 50 ) , Ft 10) 

DIMENSION FF ( 3 ) , XX ( 3 1 « ERR ( 3 ) , AA (3, 10) , WORK! IO) 

1 , Al ( 10 ) , FX(IO), TF(10) 

5 ,-FPC 10 ) , GP ( U) > “ 

REAL MASPAT 

C JET' MIXING WITH PRIMARY -AND ~SECONOA1LY~STREAMSCtF“THESAMF PERFECT-GAS 

C AND CORRECT NOZZLF EXPANSION 

X READ DATA 

ICO CONTINUE 

CALL DIFFEQ (1,X,Y,0YI 

READ! 5, 18, END = 200 ) THET A , SH AP F , V I SC , R Z E RO 

X THETA - « — INLET 8, L •'' MOMENTUM THI CKNESS * — FT 

C SHAPE = INLET B.L. SHAPE FACTOR 

C VISC - KINEMATIC VISCOSITY, C T** 2 / SEC 

C RZERO = INLET DUCT RADIUS, FT 

~ 18 FORMAT (4F10.6J 

READ! 5, 19) 

-C ' ONF CARD TO IDENTIFY THE -SOLUTION 

19 FORMAT ! 80H 

1 » 

C DUCT GECMFTRY AT INLET 

A = 3.14I6*RZER0*RZER3 

C APPROXIMATE NOZZLE AREA FT**2 

RFAO< 5, 15 ) DELTAX, XTUBE--,-TURBNO,NSUBfNGAM,XCORE, ANCZ 

NGAM=0 SUPRESS GAMMA =1 DO NOT 
NSUB ’*• I READ NFW MIXING TUBE PROFILE DATA =2 DO NOT 

15 FORMAT! 3F10.4.2 15, 2F10. 5) 

RFAD(5, 161 POO, TOO, P0l,T01,AMASSl»AMASSG, AG.RG 

PO - PS I A 

—-TO — OF GR 

AMASS - LBM/SEC 

RG FT-LPF/LBM-DEGR ; “ 

16 FORMAT ( 6F7 . 3 ) 

WRITE! 6,19) ‘ 

WRITF (6,17 J DEL TAX ,SHA°E, RZERO, XTU BE , NSUB , THETA , TURRND ,NGAH, VI SC, 

l AC.RG.POO, T00 . P01 , T01 , AMASSO , AM AS SI , A NOZ , XCORE 

17 FORMAT! ///, 5X, 7HDFL T AX = , c 10. 4 , 3 X , 6HSH APF= , F10, 6 ,4X,6HRZE R0= ,F10. 6, 

1 ?H C T , / , 5 X , 6HXTUB E = , F10. 4, 4X , 5HNSUB = , 1 5 , 9X ,6HTHE TA = ,F IO. 6.2HFT , / , 

25X,7HTURBN0=,F1C.4,4X,5HN0AM=, 15, l OX, 5HVI SC = , F 1 ... A , 9HF T* *2/ SFC , / , 

39X,3HAG=, F10.4, 26X, 3HRG=,F 10. 4, 1 5HFT-LBF/LRM-DEGR, //, 

48X,AHP00 = , F10. 4.4HPSIA, Z, RX , 4MT0C = , F10. 4 , 4HDEGR , / , 

58X,4HP01«» Fie.4,4HPS I A,7 , ax ,4HTQ1 =,F10=.4*4HOEC-R ,/ f — 

65X,7HAMASS0=,F10.4, 7HLBMZSEC , / , 5X , 7HA MASS 1 = , FI 0 . 4 , 7HL BM/ SEC ,/, 
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77 X f SH AN07*= t F 10* 5# 5HFT**2. / ,6X »6HXC0RE«,F10. A ,// ) — 

C THETA = INLET B.L * MDMFNTUM THICKNESS, FT 

C SHAPE = INLET B.L. SHAPE FACTOR 

C VISC = KINEVATIC VISCOSITY, FT**2/SEC 

C -RZFRO = INLET- OUCT R.ADIUSr FT 

C ANOZ = APPPOXI HATE N0Z7LF APFA FT**2 

C XCORF = LFNGTH OF POTENTIAL CORE/RZERO 

C DUCT GEOMETRY AT INLET 

HS c 10 

J POLY = 5 

C CALCULATION OF F(I) FROM POLYNOMIAL 

AAZ FRO = l.CJCA 

“ Al(l)= -0.0175 

A1 ( 2 ) = -8.3821 

— :A1 C 3 1 = 16.5 806' 

A1 ( A ) = - 1 2*7877 

A 1 ( 5 ) = 3.605 8 

CTR = -.05 

DC 807 1=1. NS 

CTR = CTR + 0.1 

Ftn”= Al(JPOLY) 

FP( I) = FLOAT ( JPOL Y ) *F ( I ) 

00 808~J=?f JPOLY 

JJ = J POL Y - J ♦ 1 

FU) - FU)*CTR ♦ A 1 ( J J ) 

FP(I) = FPm*CTR ♦ FLOAT ( J J )*A 1 ( J J) 

8C 8 CONTINUE - — ” 

F( I) = FIT )*CTR + AAZERO 

Tr( J7- $QRT(F< m 

8C7 CONTINUE 

WRITE (6, *03) ( F ( T ) | i M , NS ) 

AC 3 FORMAT! 5X,6HF( I) * S, 5 X , 10F 10. A , / ) 

A—*- 3. 1 A 16*RZER0*RZ ERO 

POl =P01 * l A A. 

-p 0 0 = POO* 1 A A. — 

RG =RG*32. 2 

AMASSl=AMASSl/32.2 

AMASSO=AM ASSO/32.2 

XYZ=0. 0 — 

GC TO ( 10, 11 ) ,NSUB 

10 — CALL 5UB( X YZ ,R , DR, 3,RZER0) — — 

11 JRUNGE = A 

NPART = 1 ______ — 


TO DETERMINE TOP HAT PARAMETERS 


TO DETFRMINF EFFECTIVE DUCT AREA AT END OF CORE 

CALL SUR(XCORE,R,DR, 2 , R Z E RO ) - 

THETA = THFTA ♦ X COR E*. 00 10*R ZE RO 
RCORE = R*R7FR0 - TMFTA^SHAPE 
ACORE = 3. lA16*P.CrRE*RC0RE 

AFFF ^ ACORE - ANOZ 

AM = ACORE 

A FAC = 1. 

AGG = AG/( AC-1. ) 

RHOO = POO/{ RG*TOO) - 

CORRECTION FOR INLFT DUCT PRFSSIJRE LOSS 

POC LOSS AT INLFT OR NOT 

P02f POO -SDL OSS* AM AS SO** 2. 7! 2.*RHC 0*ASC**2. ) 
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RHOO = PO2/IRC.*T00) — 

C /IRST APPROXIMATION FOP OUTER VELCCITY 

UO = AMASSO/(RHOO*AEFF I ' ----- 

TO = TOO - (AO-1. )*. S*UO*UO/( AG*RGI 

p— » — P02 + ( T O/TOOO** AOG — 

RHO=P/( RG*T01 

-c SECOND APPROXIMATION FOR OUTFP VELOCITY — 

UO «= AM A S SO / < R HO* A E c F ) 

TO = TOO - (AO-1. I*. 5*UO*UO/< AG*RGJ 

P * P02 * ( TO/TOO )**AGG 

RHO=P/( RC*T01 — 

C ASSUME PRIMARY STREAM EXPANDS TO OUTER STREAM PRESSURE 

Tt = TC1*( P/PCi J**( l./AGG) 

RH1 = P / ( RG*T 1 ) 

Ul~= SORT ( 2.*AGG*RG*( T01-T11 1 

ALIT AAO 


UJOO=trt*-UO ~~ 

T = f TO 1- TOO ) / TOO 

SOO=0. 5* l AG— 1 • ) *U JOO*UJOO / 1 AG*RG* TOO) 

C MAS = AMASS1 AM AS SO 

CPOM-* (P-POO>*AM»-AMASS1*U1 "♦ ATIASSO*UO ALIT 5 

CM = CMOM/ A *2. 

CENR-=-AHASSl+T©l >- AMAGSO*TOO 

FN=CMAS/SQRT( 2 . *CMOM* A*R HOO I 

C=UJOO*SQRT( A*PHOO*.5/CMOMI - 

W 2 * ACORE*PHOO*UJOQ 

CMAS = CMAS/W2 

CMOM = CMOM/I W2+UJ001 

CENR — = — CENR/ ( W2*TOO I 

FIRST GUESS OF PARAMETERS AFTER CONSTANT AREA TRANSITION 

-UR- = UJ/UJOO 

Y 2 ^LAMBDA 

DCR= DELTA /R “ 

UR = 1. 

Y2-UO/UJOO — “ 

PH20 =-( POO-PI*. 193 

-^TMOM= AM ASS 1*U l 

RHOC= RHO*32. 2 

WRITEI6.67) UO,RHOC,UJOO,TO,P,TMOM,Y2,PH20 ~ 

67 FORMAT( 5X .AOHCONDI TIONS AT BEGINNING OF THE TRANSITION SECTION, 

i//,9X,3HU0=,F6.0*6HFT/SHC,6X,AHRHC = ,F6i4-,-9HLRM/FT**3,6X, 

25HUJ0O=,F6.O, 6HFT/SEC, /,9X,3HT0=, F6.0 , 4HDEGR ,10 X ,2HP= , F6.0, 

39HLBM/FT ** 2 , 6X , 13HPR. I MARY MOM.= , F fc, 2 , 3HLBF , / ,5 X ,7HL AMBCA= , F6.4 , 12X 

A, 5HPH20 = , F6.2, 6HI N. H20) 

XX( 1) * UR 

XXI2) = Y2 

XX< 3 J = - SORT! SOR T! ANOZ/3. 1* 16) /R ZERO I — 

C 

DO 121 J = 1,3 — — 

121 ERRIJ) = 1. E-6*XX ( J I 

IP =1 

MAX = 120 

CAtt- VBOAAf 3,3-,-FF , X X , FB R , A A , 30£>000.,I PyMAX # 3 ,rlO»KJ 

C CHECK C ALCUl AT ION OF EN AFTER TRANSITION 
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“CMAS2 CM AS I *W2 — — 

CMDM2 * (FF(2) «■ CM0M)*W2*IJJ0G 

-EN2 = CMAS2/SQRT( ?• *CMOM2*A*ftH0Ol 

WRITE(6,6B)XX(1) » X X ( 2 ) , X X ( 3 ) , E N , F N2 

6fr-"F0RM4T(Y//,-7X,6HXX( 1 I=,F15. 10,3 X, 6MXX C 2 ) -,F15. 1 0,3 X ,6HXX <3 » = , FI 5 

1 •10,//,10X,3HEN=,FXO.A, 10X, 4HEN2=,F 10,4) 


0 * C*XX(1 ) - 

SO = $ 00 *XXm*XX(l ) 

BP-— ( AG-W )/( 2.* AG*S0 > 

UJO = XX( 1) *UJQQ 

UCENT * UJO*(l.+XX{2>> 

WRITE (6, 6G )SO,BP,UJO,UCFNT 

S5- --FORMAT! /, 10X,3HS0* , FIO. 4,11 X , 3HB P*,F lO,-4 , // ,9X ,4HUJO= »F 10, 1 » 6HFT / 

lSEC,2X,6HUCFNT=,F10.1 f 6HFT/SEC) 



Y(l) = 1. 

Y!2> *XX ( 2 ) 

Y ( 3 ) * X X ( 3 ) *R 

Y141-*! 1* Ym*Y< 2)*S0*A€/UG~H4->*Wtt/-P0O 

Y ( 5 1 = (TO l ~ TOO) /TOO 
Y (6 ) = 1,0 
Y ( 7 ) = THE T A/R ZERO 

Y!8) = SHAPE 




D=4,08*Y t 3 )* ( 1 , *Y ( 2 1 ) 

0 - XCOR E — 

X CUT ( 1) =D 

XTRO - XTUBE/R7ERO 

M 1= 2.* XTRO 

T)C 5”! = 2 , M l “ 

X0UT<1) = X0UTU-1>*DELTAX 

TF(XOUT( D.GT.XTRO) GO TO 300 

5 CONTINUE 

"3CQ CONTINUE 

M 1= 1-1 

— — — 0 C-8K 

YMIN(K) = *01 

TOt(K) = .00001 

TOL(K) = ,0001 

0 continue — 

YMI N ( 5 ) = 1. 

y^^frlT-1-- 

YYIN ( 7 ) = 1. 

YMIN(R) = 10. 

MARK ( 11 = 1 

MARK ( 2 ) = M 1 

MARK ( A) =0 

H=* 1 

CALL RUNGE ( 8 , 0, Y , TOL , YMI N , H, XOUT , MAP K ) 

GO TO 100 

200 STOP 

EN0 * 

SUBROUTINE CIFFFO (N,X,Y,nY) 

COMMON TUR BNO , CF , EN , NP AR T , CvJRU NGE ,TJG A M 

1, POOfPOl, TOOfTOl, AM ACM, AG,RG 



ALIT 

10 

-ALIT 

20 

ALIT 

30 
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— 2 - “ vAMASSl* AVASSO, T ALIT 40 

3 ,S00, CMAS, CMOM, CENR ALIT 50 

4 ,S0, BP, CM, VISC1 --- - - 

COMMON 7 A, ZB, UJO, DOR IN 

6~ , MASRAT, RZ ERO, AF AC , CFD, HD, F-, “AH “ A A 7 ErRO ,HT F F P , GP — 

7, SDL OSS, AS D 

8, AIJG l fUltUCFNT 

DIMENSION Y( 10> >BY( SOlr V(9,l), K(40I, XCUT(IOO) ALTT 

— 1 i Al(iO), TF(10) ~ 

2 , A( 70 ) » B ( 10), Ft 10) , G( 20) ,DLF ( 20) ,DLG I 20 ) 

10 ) — 

4 , DLF 1 ( 20 ) 

5 , FP( 10) ,GP( 10) t : * 


EQUIVALENCE (SO, SO) ALIT 120 

REAL HASP AT ' 

IF f N-l ) 1 * 1 , 10 ALIT 13) 

i C 0NT3 NU E - 

LM = 0 

MM=0 - ALTT 520 

MMM=0 ALTT 530 


READ! 5 , 399 ) NS — - 

399 FORMAT (312) 

~)POLY~- 5 

READ (5,401 ) ( G G( I ) , 1=1 ,NS) 

-401 ' FDRMAT( 10F5.4) — - ~ - 

READ(5,950) SrLOSS,ASn 

950 FORMAT ( 2F1 0.3) — 

WR IT E ( 6 , 405 ) 

4C5 FORMAT! IH1 ) 

HR I TE ( 6 , 400 ) NS 

4C3~ ‘ FORMAT ( 5X, 5HNS - ,1?) “■ 

4C4 FORMAT ( 5X,7FGG( I) • S , 5X , 10F 10. 4 , // ) 

HR I T E ( 6 f 4 04 ) (GCCI ) , 1 = 1, NS) 

HRITF(6,951)SntOSS, ASD 

951 FORMAT! 5'X ,~24HDUCT L OSS “C9FFF1C I 10. 3, /, t 7X, 1 2HDUCT AREA *~ f - 

1F10.3,5HFT**?,//) 

C 

DIVNS ■ 2. /FLOAT(NS) — — — - 

CIVDEL = FLOAT ( NS ) 

RETURN — 


C 

< Y1=UJ/UJZFRC — 

C Y2 = L AMDA = UO/UJ 

C Y3- BELT A/R ZERO : 

C Y4 = P /POO INITIAL 

-€ — ~Y5- ( TOCF NTf R—TOO ) /TOO 

C. Y6=T00 /TOCINITIAL 

C Y7-^~B.- L. MOMENTUM THI CKNESS /R ZF RC 

C Y8 = B. L. SHAPE FACTOR 

C" — — 

10 CONTINUE 

“4U2 'FORMAT— t4F 10.4 )— 

IF C LM» 913,913,914 

913 CCNTINUE 

0Y< 1) = 0. 

DY( 2 ) = 0. 

914 LM = 1 

f«lr/TURBN13 

Y 1=Y ( 1) 


ALIT 55) 
“ALIT 560 
ALTT 57) 
ALIT 53) 
ALIT 590 
— ALTT 600 
ALIT 610 


ALIT 620 


ALIT 640 
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Y2«Y C 2) — ' 

Y 3= Y ( 3 ) 

Y4=Y(4) 

Y5=Y { 5 ) 

Y 6=Y ( 6 ) : 

Y 7*Y < 7 ) 

ve=Ym 

Yil«YCH*Y(l) 

Y22=Y ( 2 ) *Y ( 2 ) HIT 

Y33= Y ( 3 I *Y ( 3 ) ALIT 

01VYIM./Y1 — 

DIVYll = D I V Y1 *01 VY 1 

DIVY4M. /Y4 

DIVY6 = 1. /YC 6 ) 

— $2~Y2/Il.-$0*YU+Y22) 

0S2YI=2.*S0*S2*S2*Y1*Y2 

&S2Y2 ~ 1*^( 1 *-SO*Y ll*Y22 K-^~2* *SQ*Yi^*52*S2- 

CALL $UB(X, R,DR,2,RZER0J 


C 

R=R-Y7*Y 8 * 

C NOW BRANCH TO PART 1 OR PART 2 

CO TO- Ml ♦ 20 > ♦ NRART 

11 CONTINUE 

C 

C PART 1 

IF ( Y2> 40 9tM0,411— — — 

409 WRITE ( 6» lOOO ) 

- 1000 fORM ATI 10X , 46HRECMOUL A T I ON~PR£ SENT, CALOULATIPN NOT CORRECT) 

410 OERU = 0 

- — 0YI71 * 0. — — - 

DY(8) = 0. , 

-GO TO 412 — 

411 CONTINUE 

Y7 * Y7 + .0001 
RTH = Y2*Y1*Y7/VI5C1 

CF=0.246/( f RTH)**. 266*10. *M ."678* Y8) ) 


ML IT 
ALIT 
ALIT 
ALIT 


412 


C 

C 

C 

G 

C 

C 

C 

c 

13 


CFD = CF 
CCNTINUE ' — 
DIVY3=l./Y3 
E~ g E 


Ml. ♦ 1,5*11# - 2. 718**1-1. 1*Y2) ) 1 


HO 


*1 
* 2 
*- 3 - 
= 4 
= 5 
7 

- $ 


POO CONSTANT — 

MOMENTUM INTEGRAL 

“CONTINUITY-INTEGRAL 

ENERGY INTEGRAL 

MOMENTUM HALF INTEGRAL 

B.L. MOMENTUM EQUATION EQN 
B.L. MQMENT-OFMOMENTUM EQUATION ” 


CONTINUE 
DCNR=Y33/( R*R) 
CALCULATION OF 

2A = 0 . 

7 B = 0. 


EN 


PART 1 


CTR .05 

DC 369 I = 1 1 NS 
"GTR * CTR" 1 — “ 
Y2FI = Y2 + F { I ) 


ALIT 

ALIT 

ALIT 

ALIT 

ALIT 

ALIT 


ALIT 

ALIT 


6 5 J 
64 > 


Ml 
710 
4 70 
7 2) 


7 31 
7'0 
7 5 J 
763 
773 
7 80 


7 93 

3)0 
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O (-* 


FS = TF( I) - - - - 

ANA = Y2FI 

ANB - Y2FI*Y2FT "" 

0 * 1# ♦ Y(5)*FS -S0*Y11+Y2FI*Y2FI 

ZA s ZA + D1VNS *ANA *01 VD 

1 *CTR - 

ZB * ZB ♦ DI VNS*ANB*DI YD 

1 *CTR - 

369 CONTINUE 

S2 - Y2/( 1.-S0*Y1I*Y2*Y2) — 

AM AS = ZA+OONR ♦ (l.-D0NP)*S2 

AMOM ZB*DONR + ( 1 . -DOUR ) * Y2*$ 2 " * - 

FN = AMAS/SORT( AMOM «■ BP* ( 1 . - 1. /Y ( 4)) *01 VYU ) 

EN * EN* # 707 - 

RLOCAl = R*PZERO 

RLOCAL ■*’ R LOG AL*R LOCAL " 

MASRAT = A M AS*Y4*POO*Y1 *U JO* RLOC A L *3 . 14 1 6 / ( RC* TOO* AMA SS 1) - t 


W<l f n=Y22*Yl 

W(i, 2)=YU*Y2 

W( 1 , 3 > O. 

W(1,4)=BP*( 1.-$0*Y1 l*Y22l*DIVY4 

— * W ( 1 1 5 ) = 0. 

WU f 7l = 0. 

W(l,3> -* 0. 

VCltl) = 0. 

J * 1 

1\ J = J+1 

71 *-0 

7 2=0# 

Z 3=0* — 

Z4=0. 

z 5-$.- 

Z6=0. 

CTR = -.-05— 

DO 66 I=l»NS 

CTR = CTR ♦ .1 

FS = TFIU 
F I = F ( I) 

Y 2F l = Y2 + F ( I ) 

J J * J- 1 

GO TO ( 270 f 70 1 67 ) ,JJ 

2 70 CONTINUE 

AN=Y2F I *Y2 F I 

D2N=Y2Fl^2. 

05N=0 

GO TO 68 

70 AN= Y2F I 

D2N=1. 

D5N=0. . 

GO-TO-63— 

67 AN=Y 2 F I * ( 1 . +Y5*FS ) 

D2N= l • *Y5* FS — 

D5N=Y2FI*FS 

*0— 0 = 1. +Y5*FS-S0*Y11*Y 2F I* Y2F I 

0 IVD= 1 * /O 

— TTIVDD-TTI VD*DI VD 

D2D=-2.*S0*Y11*Y2FI 


ALIT 8 V) 
ALIT 853 


ALIT 863 
ALIT n ) 
ALIT 923 
ALIT 933 
ALIT 9 ^ j 
ALIT 95) 
ALIT 963 
ALIT 97) 
ALIT 93) 


AL IT 1 ) 1) 
AL I T 1 3 2 3 
ALIT103) 


ALIT107) 


” - ALIT! UU 
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01D=-2.*SO*Y1*Y2F I*Y2FI 

05D=ES 

Z 1 = 7 1 +DI VNS*/'N*DI VD 

1 *CTR 

2 2 = 7. 2 -n 1 VNS*AN*D1 D*D I VtH> 

1 *CTR 

Z3=Z3+DIVNS*D2N*OiVD 

1 *CTR 

£4-2 4-0 1 VNS* AN*B20*0 1 VDO 

1 *CTR 

25=Z5*DI VN$*D5N*D I VO — 


1 *CTR 

-Z6=Z6-Ol VNS*D50*0TVO0*AN 


1 *CTR 

66 CONTINUE 

IF ( J-2) 210,211,210 


1LTT1250 

SUT12M 


211 


-W(j!n=OONR*i2.*21*OIVYl>-721Ml .*DONR »*12.*Y2*S2*DIVYl*Y2*OS2Yl I 

W IJ « 2 1 = DON R* { 7 3 +Z 4) +(1.-D0NR)*( S2 + Y2* CS2Y2 I ^ _ 

H (J, 3 1= 2.*0CNR*( Z 1- Y2*S2 )*D! VY3 “ 

WU,4» = ID0NR*I21-Y2*S2H-Y2*S2»*DIVY4 _ 

-1 + 0P*Of VY4*DI VY1 1 ' ‘ 

W(J. 5)=DONR*IZ5*Z6) 


HU, 7) — 2.*S2 + Y8/R*Y2 — " 

H(J, 9) =-2.*S2*Y7/R*Y2 

V(J,ll=-2.*OR/R+Y2*S2 

-210 HI Jil )=0CNR<'<Z1*DIVY1 + Z2) + U.-00NP»*( S2*0I VYU0S2 Y1 1 
HI J,2) =OONR*< Z3+Z4) + ( l.-OONR ) *DS 2 Y2 
HU » 3 1 = 2 . ♦DONfi'M 71 _ S2)*DIVY3 


H(JiA) = (D0NP 1 *(Zl-S2>+S2M i 0IVY4 

HIJ, 5 >=D0NR*(Z5+Z6) 

W(J, 7) =-2.*S2*Y8/R 
HU, 3) =-2.*S2*Y7/R 


VU,1)=“2»*DR/R*S2 

212 ^-CONTINUE 

IF ( J— 4 ) 71, 6°, 69 

'69 — CONTINUE 

J = 5 

-226 1 = 0 

21 = 0 

2-2 - 0 

23 = 0 

Z4 = 0 

25 = 0 

— — : — 26 = 0 • 

27 = 0. 

28 = 0. 

Z 9 = 0. 

Z 10 = 0. — 

711 = 0. 

Q1 = C . - 

02 = 0. 

Q3 = 0. 

Q4 = 0. 

HI = 0. 

R 2 = 0. 

R34 = 0. 

R56 = 0. 

R79 = 0.' 

RIO = 0. 


ALITl 

4 90 

ALIU 

53 ) 

ALIU 

51 > 

M. I T 1 

530 

ALIT 

19 

ALIT 

29 

ALIT 

— 3 > 

ALIT 

A) 

ALIT 

50 

MIT 

5 3 

ALIT 

70 

ALIT 

80 

ALIT 

99 

ALIT 

no 

ALIT 

II ) 

ALTT 

120 

ALIT 

131 

ALIT 

143 

AL I T 

159 

ALIT 

1 63 

ALIT 

173 

ALIT 

130 

ALIT 

190 

ALIT 

2 )0 

ALIT 

219 
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rii * o. 

vi = o. 

v? = o. — 

V3 = 0. 

V4 * — 

V5 = 0, 

V6 = 0. - 

V7 = 0. 

V9 = 0. 

VI 0 * 0. 

vn -«--©• 

CTR = -• 25*DIVNS 

224 1 *1 + 1 

CTR = CTR ♦ 0# 5*0 IVNS 

FS = TF< I ) 

Y2FI * Y2 4 F( I) 

FDER— FfM M — — — 

221 AN “ Y2FI 

02N= l. 

D5N=0. 

223 D = l. 4Y5*FS-S0*YU*Y2FI*Y2FI — 

DIVD=1./D 

0 IVDD-0TVD*0I Vf3 

D2D=-2.*S0*Y11*Y2FI 

01D=-2.*S0*Y1*Y2F I + Y2FI 

D5D=FS 

CCC *"DIVNS*CTR*.5 

0Z1 = CCC*AN*DIVD 

0 Z 2- =^CCC * A N * D 1 D*0 T VDt> — — 

0Z3 = CCC*D2N*DIVD 

0Z4 - -CCC*AN*D2D*DIVD0 — 

1)75 = CCC*D5N*DIVD 

0Z 6 - -CCC *AN*05D*DTVDD 

D710 * FDFR*CTR*CCC*DIVD 

— DZll-«-CCC+C.TR*Y2P! *01 VDD* 2.* S0*Ytl*Y2F I *FDFR 

VI ^ Z1 ♦ 0.5*071 

V? * Z2 ♦ 0.5*DZ2 

V3 = Z3 4 0. 5*DZ3 

V4 = Z4 4 0.5*DZ4 

V5 = Z5 4 0. 5*0 Z5 

V6~ Z6" 4 0. 5* 07 6 — — 

V10 = Z10 4 0.5*0710 

Vll -- Z1 1-4 .5*02 11 — 

Z1 = Zl 4 DZ1 

12 ^ Z2 4 DZ2 — 

Z3 = Z3 4 DZ 3 

t*~=- Z 4~ 4—02 4 

Z 5 * 75 4 OZ 5 

-Z6 Z6 4 026 

710 = Z10 4 DZ 10 

7 11 = 711 4 om - 

492 CCC = 0 1 VN S *C TR*C TR 
1 *.5 

ATAU * .37 7 

BTAU * .0333 — - — 

Q1 = Q1 + CCC*Y2FI*DTVD 

1 *Y2FT 

Q2 s 02 ♦ CCC*Y2FI*0IVD 

Q 3 ~03 ~*~CC C * C TR *Y2FI*FOER*T)IVO 

CCC = -CCC /CTR 


AL 

IT 

2 7 ) 

AL 

IT 

2 3) 

AL 

IT 

740 

AL 

IT 

75 ) 

AL 

IT 

750 

AL 

IT 

770 

AL 

IT 

790 

AL 

IT 

7 90 

AL 

IT 

3)'< 

AL 

IT 

31 J 

AL 

IT 

370 

AL 

IT 

3 30 

AL 

IT 

340 

AL 

IT 

350 

AL 

IT 

3 39 

AL 

IT 

39 3 

AL 

IT 

43 ) 

AL 

IT 

470 

AL 

IT 

4 3 ) 

*L 

IT 

440 

AL 

IT 

4 50 

AL 

IT 

440 

AL 

IT 

4 7 j 

AL 

IT 

490 

AL 

IT 

53 3 

M 

IT 

-510 

AL 

IT 

570 

AL 

IT 

530 

AL 

TT 

540 

AL 

IT 

55 3 

AL 

IT 

539 

AL 

TT 

5 9 9 

AL 

I T 

V)y 

AL 

IT 

sk 

AL 

IT 

5 70 

AL 

IT 

633 

AL 

IT 

64* 

AL 

IT 

670 

AL 

IT 

630 

AL 

IT 

690 

AL 

IT 

733 

— m 

IT 

713 

AL 

IT 

7 70 

AL 

IT 

7 30 

AL 

IT 

769 

AL 

IT 

790 


ALIT 303 
4LIT 310 
ALIT 370 
ALIT 300 
ALIT 9J0 
ALIT 910 
ALIT "9 7 ) 
ALIT 94J 
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R1 = R1 CCC*V1*FDER 

R2 = R2 *■ CCC*V2* FDER 

fm = R34 ♦ CCC*I V3 ♦ V4)*F0ER 

R56 = R56 + CCOMV5 + V6)*FDER 

RIO * - R10~“~CCC*V1.0*FOER “ 

R 1 1 = R1 l-CCC*V 1 1 * F DER 

IF- I I-MS 1 224,225,2 25 

225 CONTINUE 

VMJ.1I = ( Q1 ♦ Rll*OIVYl ♦-«?- 

WIJ.21 = 02 ♦ R34 


ALJT 95-) 

ALIT 950 

ALIT 970 

ALIT 990 

— "AVI T 1-3 0? 

ALITIOIO 

ALIT1J2Q 


9M-J, 3> t 03' ♦ "RlO '^-Rt th*OtVY3 ” 

y(j,4> = R 1 *DI VY4 + 8P*DIVY4*DIVY11*.3333 



W1 J,5 1 = R56 

WIJ.71 = 0. 


■ 

— ~W I J , 8 1 = 0.- 

VlJ.l) = E*DI VY3* AT AU 


’ 

1 


245 

continue 


246 

—“CONTINUE 


247 

CONTINUE 



947 CONTINUE 

C B.L. MOMENTUM EQUATION 


Vft7,i)-= i^.-r-Y817Y7*0IVYl*Y2 

HI7.21 = I 2. ♦ Y&1*Y7 

Vfl 7r3 1 = 0. " 

W 17, 4 1 = 0. 

kL7 , 5 >~ *"©• 

HI7.61 = 0. 

W(7r7»- * “Y2 


W(7,8) = 0. 

VVTri l 5*CF*Y2 “ 

C B.L. MO MF NT- OF -MOMENTUM EQUATION 

HI8r21 ~.5*Y8*IY8 ♦ 1. 1*1 Y8*Y8-1 

W(8,1J * W ( 8,2 >*DIVY1*Y2 

VMS, 3 V 0 '. ” 

W 1 8, 41 = 0.- 

— WI8.51 = 0. 

W(8,6) = 0. 

il ( 8 f 7 1 ~Om 

W(8 81 * YZ 

V( 8* 1> 1 Y8*Y8 — l.-)*Y2/Y7*( .5*Y8 , t , CF-.06*( Y8-1. ) / ( ( Yt)+3. 1 

1 *RTH**.lll 

I F- I YI21) 20<h 201* 201 — 

200 W(l,l 1 * 0. 

HI 1* 21 = 0. 

WI7.1) = 0. 

M 17,21=— 0. ’ ” 

V ( 7 , l ) = 0. 

HI 8,1)-- 0 . — — 

W 1 8 1 2 I =0. 

VI 8,11 = 0. 

C0NT1NUF 

-I-F— IY8-2. 1 46?T*61f~581 — 


2C1 


-AL IV 870 
ALIT 883 
ALIT fl90 


ALIT 9 JO 


481 H ( 7 , 1 1 = 0. 

-HI 7, 21 -= 0. 

W 1 7, 7 1 = 1. 
VI7.11 


HIS, 11 - 0. 
HI 8, 21 = O. 
HIS, 81 = 1. 
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V(8,l> = 0. — 

- 0. 

V(l,2l * 0. — 

482 CONTINUE 

C COLL APS E~ FROM 8X 8 TO TXT MAT R I X " 

DO 687 1=1,8 

w {6? 1 1 = w(7,n - - 

W f 7* I I = W( 8, I ) 

687 CONTINUE 

V(6,ll = V ( 7 » 1 ) 

w,n r- V (t f n 

c 

DC 688 J = 1 , 7 

MUf 6} , VMJ,7) -■*■■■ 

W ( J, 7 ) = W C J « 8 ) 

*88 CONTI NUE — ? 

C 

C ALIT2153 

C SOLVING SIMULTANEOUS EQUATIONS 8Y SIMQ SUBROUTINE — ~ — ” ALIT2160 

C ALIT2170 

NN^= 7 — 


DO 101 J = 1 , NN 

B ( J ) =V{ J , 1 I 

DO 101 1*1, NN 

U * 1 ♦ NN*t J-l) — 

A 1 1 J ) =W ( I, J) 

— iO I ~ CONTI NU E- 

CALL SIMQ (AtBfNNfKS) 

DO 18 1 = 1, NN — 

IS DYIl)=B(I) 

DV( 8) = DY ( 7 ) 

0Y(7> = DY ( 6 ) 

DYt 6 >-0* — — 

249 CONTINUE 

JRUNSE= JRUNGE + 1 

C CHECK FOR EQUATION CHANGE 

IF ( JRUNGE-5 I 14,15,15 

15 JRUNGF=0 

If IYL3T- R1 — 14 r 1 6,1 6 — 

14 RETURN 

16 — X0UT<]OO)«X — — 

CALL PRINT (N, XOUT, Y,DY , 100) 

Y<31=0. 

Y3 = VC 3 ) 

NPART*2 — 

BPPT1 = BP 

S0PT1 = SO 

DLSTAR = Y8+Y7 

DELTA - 0LSTAR*(Y8*l. >/( Y6-1.1 - 
POW = 0.5*(Y8-1.) 

' R = R + OL STAR— — 

C EVALUATE NEW VELOCITY PROFILES 

CTR = - # 05 

DO 876 1=1, NS 

— — CTR ■» CTR ♦ 0.1 

IF Il.-CTR-DLSTAR /R ) 875,875,874 

-875- FH I— 0. 

TF( 1 1 = 0. 


ALIT2170 
4LIT2230 
AL IT2 2 10 

ALIT2270 

ALIT2230 

-AtlT2 24-y 

ALIT??5.J 

* — ' ALIT2 261 

4LIT227 j 


ALIF7230 

AL I T 21 40 

ALIT2340 
AL 1^2 3 50 
MIT2363 

ALIT2370 

ALIT7333 
A L IT? 370 
UIT24JQ 
ALIT2410 

ALIT24M 

ALIT243 ) 
1L I T 2440 
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874 


873 


840 


861 

860 


876 


CO TO 840 - 

CTRl = CTR / 1 1.-0L ST AR/R ( 

8 C 1 1 = A 1 ( JPOLV ) 

00 873 J=2 t JPOLY 

JJ = JPOLY - J ♦ 1 

Fill = F ( 1 1 *CTR 1 * All JJ^) 

CONTINUE 

F( II = F( I 1 *f,TR 1 ♦ AAZERO 

TFIII = SQPT ( F ( 1)1 

BL= ( R/DELTAM l.-CTR) )**POW 

— |F ( BL-IV)~860t 860 f 861~ 

BL = 1. 

CONTINUE 

Fill = (F( I) ♦ Y2*BL)/(1. ♦ 
-cm GGIII - Fill 
CONTINUE 

HR I TE I 6 * 40 3 1 (Fill, I-l.NS) 


Y2) 


4C3 

406 


FORMAT t 5X,7HF2C II 1 St IX 1 10F 10.4 t //) 

WRIT E ( 6 i 40 6 ) ( C( T ) t I =1 >NS) 


407 


-877- 


FORNATI 5Xt 7HG2 ( I ) ’ St 5X , 1QF 10.4 , //I 

FORMATl"^ 07 * //4X.7HX/RZER0 , 5X , 4HARE A , 5X , 4HPH2C 

8/ TOO^UC ENT t FT/ SEC I -- UR 7HAUGMENT, 6X, 5HGANMA 

3NTCOEGR) TGWALUDEGm ) 

-t- ss "NS"*-l 

00 877 I -2 tt _ , t4 

~fP( I ) -.5 *OIVDEt*I Ff !♦! t— Ft 1-1 > 1 
GP( I ) =.5 *DIVDEL*( G(I + 1 ) -G( I-t )> 

“CONTINUE 

FPfl) = OtV&EL*fF(2 ) - FUJI 
” GPt II s DT VO EL* (0(2) - Of 1 M 
FP(NS) - D I VDFL*( UNS) - FfNS-1)) 

CP(NS) * CIVDEL*<C(NS> 0(NS-l 11 
Y 1 = Yl*( 1 . «■ Y2) 

-Yll - - Yl + Yi ’ 

VC1I * YIDN1. + Y2) 

- OtVYll = 1./V1.1 

E2I = E / ( 1 • ♦ Y2) 

— E21 = ?.*E 2 I 

CFOI = C F / ( ( 1 * ♦ Y21*(N ♦ Y2>) 

0 # — 


f 5 X t 3 
,3X, 


5HTOWALL 

25HTOCF 


“Y2~*- 

Y ( 2 ) 
Y 7 - 

Y ( 7 ) 

Y 8 = 

Y ( 8 ) 


= O* 

0 . — 
= 0 . 
0 . — 
= 0 . 


0 — 

c 

-c— 

c 

c 

- 20 - 


PR 


C 

--C- 


PART 2 

TURRULFNT PRANOTL NUMBER 
PR ■* 1. 

yb^gahma 

— CONTINUE 

Y23=Yf 2 >*Y< 3) 

SO = S0PTI/YI6I 

0S0Y6 = -S0/YC6) 

— ftp « BPPT1*Y(6) 

PH AS 6— OUT— T WE -6 *t~r “M S°L AC EH ENT— THICKNESS -I NFER I TED FROM PART 

0YI7J = 0. 


—At IT >48) 
ALIT2463 
-At IT 247 ) 
At I T 248 3 

At I T 2 5 n 

At 1 T2320 
ALt T 2530 
ALIT2543 
At IT 23 5 J 
At IT 2360 
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oortnrtnn 


ov(a» = o. 

C CALCULATION OF EN PART 2 

ZA = 0* 

ZB = 0. 

-CTR » -«05 . 

00 370 I = 1 , NS 

CTR --CTR + 

YFGG = Y2 «- F( I ) + Y3*GCI) 

FS - -TFU> 

ANA = YFCG 

AN 8 = YFGG*YFGG 

0 = I. ♦ Y5*FS -SO*Yll*ycGG*YFGG 

01 VO = 

ZA = ZA + DIVNS*ANA*DIVD 

I 

ZB = ZB + OIVNS*ANB*OIVD 

t~ *CTR _ 

370 CONTINUE 

a has = za — 

A MON = ZB 

= AMAS/SOPTI AMOM ♦ B P* ( 1. - 1, /Y ( At ) *01 VY11 ) 

EN x EN*. 707 

R LOG Al— *~R-*R ZERO — — 

RIOCAL = RLOCAL*PLOCAL 

MASRAT = ANAS* Y4*PQ0*Y1 *U JO*RLOCAL *3. 1416/ I RG*TOO*Y6*AMASS1 )- 1 
UM = AM AS* Y1 

RW = UM*2.*P/VISC1 

CFOF = AMA$*AMAS*.048*RM**(-. 20 ) 

_____ _ 

CFO » CFDF 


J=I CONTINUITY 

J=2 - ENERGY 

J=3 MOMENTUM 

^*=4 " MOMENT OF- MOMENTUM ECUAT ION 

J=5 T* /T = UJ'/UJ 

J=6 DILAMBDAI/OX = 0. 

J=0 

~86 J=J*1 

87 Z1=0 

7 2-0 — — 

13=0 



Z 5=0 

Z 6=0 

Z 7=0 



Z 9 = 0« 

CTR = 

1=0 

82 

CTR = CTR + .1 

F S -=- T F ( I ) 

YFGG=Y2 *F( I )+Y3*G( I ) 

GO TO ( 81, 83 t 41 1, J 

81 AN=YFGG 

D2N=i. 

D5N=0 

03N=Gm 

GO TO 85 


ALIT257 j 
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-83 AN^YFCG*! 1. *Y5*FS 1 

D2N=1.*YS*FS 

D5N=YFGG*FS — ~ 

03N=G( I)*( 1 • +Y5*F S ) 

GO TO 85 " 

41 AN=YFGG*YFGG 

- 02 N= 2 .*YFGG 

03N=2.*YFGG*G1 I) 

05N=0 

IF ( Fill > 720,720,721 

7 20 AN-=~-AN 

02N = -02N 

03N = -03N 

721 CONTINUE 

— 85 - 0=1. *Y5*FS-S0*Y11*YFGG*YFGG ~ ’ 

DIV0=1. /D 
0IV0D=01V0*DIVD 

01D=-S0*Y 1 *2. * YFGG* Y FGG 

020=-SO*Y1 1*2. *YFGG 

030 =' -S0*Y11*2.*YFGG*GI II 
— 050=FS 

06D=-0S0Y6*Y1 1*YFGG*YFGG 

27=2VrOI VNS*AN*01 VD 

1 *CTR 

— — "72*22-01 VNS*AN*01 0*0 I VDO 

1 *CTR 

28=r3*0IVNS*D2N*0IVD 

1 *CTR 

I VNS*02D*0TVDD* AN " 

1 *CTR 

Z5=Z5*0IVNS*D5N*0IV0 

1 *CTR __ 

28=76-0 1 VN$*AN*D5D*01VDD 

1 *CTR 

2-7=27*01 VNS*D3N*0I VO 

1 *CTR 

78=78-01 VNS*AN*060*0IV0O-* 

1 *C TR 

l<) 29 - -OIVNS*0I VDD*AN*D30- 

1 *CTR 

IFH-NS 182,601,601 

601 WIJ, 11=21*01 VY l +Z 2 

W< J* 21=Z3+Z4 

W ( J , 3 1 = Z7 + Z9 

W(J,4l=Zl*D1VY4 

W(J,5)=Z5*76 

WU, 61 = 78 

IF (J-2J 602,603,602 
-602 WIJ.61 = W I J « 6 1 - Zl*DIVYfe 
603 CONTINUE 

VIJ, 1 J=-2. *DP/R*Zl 

IF ( J- 3 1 8 6, 70*3 ,700 

-7(9 CONTINUE 

WI3.11 W ( 3, 1 1 * Z 1*01 VY1 

AUG1 = Z 1 

W ( 3 ,4 1 = WI3.41 * 8 P*D I VY 11 *DI VY4 

V ( 3 , 1 1 = V ( 3 , 1 1 - CFD/R “ 

J * 4 

1 ~ 

Z1 = 0 
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72 * 0 

Z3 * 0 

Z4 - O 

Z5 = 0 

Z6-«~0 — 

11 = 0. 

ze = o« 

7 9 = 0* 

— Z 10 * 0. 

ZU = 0* 

— 01 « o. 

02 = 0 • 

03 - 0, 

Q4 = 0 • 

R 1 « - O* — — 

R2 * 0. 

-R3V ^-0* — 

R 56 = 0. 

R79 * 0* — 

RIO * 0* 

R 1 1 * O. 

RTAU = 0. 

VI - 0. 

V2 = 0. 

V3 * 0. — 

V4 = 0. 

V5 = Om 

V6 = 0# 

V7 ^ O — 

V9 = 0. 

VI J = o, 

VII * 0, 

E2FMAX = -.5*CFO**45/(FP(5> + GP<5H 

HC = ( 1 • -Y33 ) *E2 I + Y33* E2FMAX 

CTR=-.25*OIVWS 

924 1=1+1 

t;TR = CTR ♦ 0« 5*0 TVNS — 

FS = TFU) 

Y2FI * y2 + F( I) 

1 ♦ Y3*G( I > 

FOER - ~FPM ) + Y3*GPft> 

GGOER = FP( I ) + GP( I ) 

— — — E2F S -.5*CFD*CTR/GGD£fl 

E = fl.-Y33)*E2t + Y33*E2F 

AN = Y2FI - 

D2N= 1 « 

05N=0« — — 

D6N = G{ I ) 

D*l« +Y5*F$-S0*Y11*Y2FI*Y2FI 

DIVf)= 1 • / D 

O T VD0=DIVD*CIV0 

D20=-2#*$C*Y11*Y2FI 

Dl0=-2« *50*Y1*Y2F T + Y2F 1~ 

D5D=F$ 

D60 = -SO*YU*2.*Y2FI*Gm 

CCC = 01 VN $*CTR* # 5 

DZ1 = CCC* AN*D IVD 

072 = -CCC* AN* DID* D IVDD 

DZ3 -“CCC*02N*DIVD 

DZ4 = -CCC*AN*D2D*DIVD0 


- MIT 3 3 
ALIT 60 

At- IT V ) 

ALIT 33 

^1!T -90 

alit m 

— ALTT 1 1 i 

ALIT 1 2 i 

4LIT ly} 

— ALIT 140 

ALIT 1 30 

At IT-160 

ALIT 170 

* 4L-IT-130 

ALIT 190 

AtFF-2 >0 

ALIT 210 

* L! T 22 , } 

ALIT 230 


ALIT 240 
ALIT 230 

'ALT T 2 60 

ALIT 270 
"ALIT ?30 
ALIT 290 

ALT T “3)0 

ALIT 31) 
“ * ALIT 320 


ALIT330 

- - - ALIT 350 

ALIT 3 6 j 
ALIT 370 


ALIT 410 

At IT “420 

ALIT 430 

ALIT 45) 

- ALIT 460 

ALIT 470 

ALIf 430 

ALIT 490 

-“ALIT 5)0 

ALIT 510 

—ALIT 520 

ALIT 530 

ALIT 540 

ALIT 550 
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075 = CCC*O5N*DIV0 “ 

076 = -CCC * AN* 050 *0 I VO D 

- 077 = CCC*D6N*niVD 

0Z8 - -CCC*D.IVD0*AN*D6D 

02^- 3r -~CCC*D6C*AN*0 I VOO “ — ~ 

0710 = FO£R*CTR*CCC*niVD 

OZll = CCC*CTR*Y2FI*DIVOD*2.*SO*Yll*Y2FI*FDER 

RTAU = RTAU-E*FOER*CCC*DIVD 

VI = 71 ♦ 0.5*021 

V 2 * Z 2 ♦ 0.5*072 

V3-*-Z3 t 0.5*073 ” 

V4 = Z4 + 0. 5* DZ4 

V5 - 7. 5 0.5*07.5 

V6 = 76 + 0. 5*026 

— — V 7 w 2 7 ♦ 0.5*027 

V9 * 29 + 0.5*029 

Vifr-*“MO * 0.5*DMO * 

VII = 211 * .5*0711 

Zi -=-21 * 021 

72=22*07 2 

— 1 3 -= _ 23 + OZ 3 — “ 

24 = Z4 + DZ4 

— 25-^-Z5~* 07 5 ' 

76 = Z6 * DZ6 

2 7 * 77 * 077 — 

Z8 = 28 + 028 

Z 9 x 79 + C2 9 

710 =710 * 0710 

211 = 211 + 0711 

CCC = DIVN$*CTR*CTR 

1 *.5 

01 * 01 * CCC* Y2 F I* 0 1 VD 

x +Y2FI * 

Q2 * 02 * CCC*Y 7 F I* D I VQ 

— Q3-^Q3-^-“CCC*CTfr*V 2f I * F OCR *0! VO ~ 

Q4 ■ Q4 * CCC*Y2F1*G( IJ*DIVO 

CCC -*’-CCC/CTR *" 

R 1 = R 1 * C CC* V 1* FDER 

R2 *" R2 + CCC* V2* FOER ' 

R34 * R34 * CCC*(V3 * V4t*FDER 

— r^R56 — *'~CCC* ( V5 " > ^V6 ) * FDER — 

R79 = R 79 * CCC* ( V7 * V9)*F0ER 

RIO — = RlO~ to “CCC*V10*FOER — “ 

R 1 1 = Rll-CCC*V11*FDER 

IF (I-NS) 924,925,925 

925 CONTINUE 

W t trt) -*-t Ql -^Rl W»TVYi~*-« ~ “ 

W(J«2) = 02 ♦ R 34 

W(J,3>-= R79 ♦-Q4 

W(J,4! = R1*0IVY4 + BP* D I VY4*0! VY 1 1*. 3333 

Wl J ♦ 5 ) = R 56 — — 

W ( J , 6 ) * 0. 

VtJf-n-— -RTAU’/ft 

1 - { 03 * RIO ♦ Rin*OR/R 

2 — -CFO*. 5 “ 

W{5,1} = DIVYl 

VM 5 , 2 ) = 0. — 

VMS, 3) = 0. 

W<5r4l-^ ~Oi— 

W ( 5 , 5 ) = -1./Y5 


ALIT 56) 
ALIT 573 
AL I T 5 30 

~M\T 5 10 
ALIT 630 


MIT 61) 
ALIT 63) 
— ALIT 633 
ALIT 64) 
ALIT 653 
ALIT 663 
ALIT 670 
ALIT 630 
ALIT-690 

ALIT 7)0 

ALIT 710 
ALIT 723 
ALIT 21) 

At IT ~7 40 

ALIT -750 
—ALIT 760 

-ALIT- 773 
ALIT 730 


ALIT 9.70 
ALT* 919 

- ALIT 920 
ALIT 91) 

At IT -940 

ALIT 950 

- ALIT 963 
ALIT 970 
ALIT 93) 
ALIT 990 

ALTT103) 

*4L IT 1)13 
AL I T 13 70 


Atm-^so 
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W<5,6>= - 1 . /Y6 * 

V(5,l) = 0. 

31 W (6 , 1 ) *“0 « 

W { 6t 2 ) * 1. 

W(6,3I~- 0^ 

W(6,4) = 0. 

W ( 6, 5 > ^ 0.- 

W (6,6) * C. 

v ( 6 # u--»- 0# - — — 

30 CONTINUE ALIT4160 


168 

169 

MM=MM+1 


A L I T 20 5 1 

CUNT INUE 
CONTINUE 



-C 

c 

~r 

SOLVING SIMULTANEOUS EQUATIONS BY SIMQ SUBROUTINE 




m = 6 

on 102 J«1 , NN 

B(J)=V(J,1 ) 

DC 102 1*1, NN 

I J = I ♦ NN*t J-l ) 

A( I JJ-Wf 1 T -J ) 

102 CONTINUE 

CALL 5 IMQ (A,B,NN f KSI — 

DO 103 1*1, NN 

— 10 3 OY(II«B<U — — 

165 CONTINUE ALIT232C 

* 16 6 C ON T I NU 6 A L-I T 23 30 

RETURN 

END = — 

SUBROUTINE PRINT ( N , XOUT , YOUT ,D Y , J ) 

DIMENSION YOUT ( 10) ,X0UTC200J t f)Y < 501 r*< 10 ) , YS AVE < 10 > 

1 ,F( 10) , A 1 C 10) , TF( 10) 

5 ,FP( 10) , GPC 10) 

COMMON TUR PNO , C F , EN , NPA R T f C , J RUNG E , NG AM 

IV POO , P01 , TOO , TO 1 , AMACH , AG ,RG 

2 , AM AS S 1 , AMASSO, T 

3 vSOO, CMAS, CMOM, CENR — — — 

4 ,S0, BP, CM, visa 

CCMMON-ZA, ZB,-UJO,- DORTN 

6 ,MASRAT,PZFRO,AFAC, CFD , HD, F , Ai, AAZERO , TF , FP, GP 

7 , SDL OSS , A$C 

8 , AUGltUl.UCENT 

REAL MASRAT ■ 

X=XOUT ( J ) 

C ALL SUB ( X , R , DR ,2 ,RZER0) 

C CHECK FOR INITIAL PRINT 

~ 2 IF « J-l) 5,5,100 

5 WRITE (6,50) 

C PRINT FOR ALL PARTS 

lOO UJ-YOUT ( 1 ) *UJO 

PH2Q = Y0UT(4) *P0O*. 193 

! POO*. 193 — 

U0= YOUT ( 1 ) *YOUT( 2 ) *UJO 

SUM= YOUT ( 1 ) *U J 0 + U0 

UR* SUM /UC ENT 

T0WALL*Y0UT(61*TQ0 * 

TOC ENT*TC< WALL* YOUT ( 5)+T0WALL 
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AREA =R**2* 

GO TO (29,30>,NPART 

2<? WR!TE(6,55)X, AREA, PH20 , t)0 , SUM, UR , 


YOUTt 2 ) , YOUTl 3 ) * TO CENT f TOWML » 


1 Y OUT (7) ,YOUT<8) 

• GO TO 31 ’ “ 

30 AUG2= 3.1416*1 <YOUT(A)*POO)/BP)* (YQUT(n**2.>*AUCl 

l*((R*RZER0)**2. ) ~ 

AUG=(AUG2/( AMASS 1*U1) ) 

*RITE<6 f 56>X,ARFA, PH20,YQUT<6) , S UM, UR , AUG, YOUT ( 3 ) , TOCENT .TOWAL L 


31 

A 1 
-AO 
50 


CONTINUE 

TF^f 0-100 V AO, AT , AO 
WRITE (6,60) UJ 
RETURN 


IJCE 

TOW 


— 55- 


56 


FORM AT ( / / / , AX , 7MX/RZ ER 0, 5X , AHA REA, 5X ,AHPH20» 5 X ,2AHU0 (F T /S EC ) 

1NT{ FT/S EC ) , AX , 2HUR, 5 X , 6 ML AM9DA, 6X , 7HDELTA/R , 2X , 3 4HTCCE NT ( DE GR ) 
2ALL(DEGR) THET A/RO , AX , 5HSHAPE ) 

FORMAT ( 3X , F7*3,4 X , F6* 3 , AX ,F 6. 2 ,5 X,F7^ T~,-5^< ,F7* I t AX , F7. A , AX , F7 3 3 , 5X 

1, F7. A,5X, E7. 1, 5X, FT, l r 6X,F7. A,AX, F5.2 ) 

— FORMAT ( 3X,F7.3,AX,F6.3, AX , F 6* 2 , 5 X , F7. A , 5X , F7. 1 , AX , F7. 4 , AX , F7* 3 ,5X 


1, F7.4,5X, F7.1,5X f F7. 1) 

-FORMAT ( // 52H OELTA/R ~ 1 --DIFFERENTIAL 

11 5X , 5HCL = ,F10.5//> 

END ~ 

SUBROUTINE calxfg U,X,FF,GG) 

COMMON TUR BNO, CF, FN ,NPART ,C , JRUNGE ,NG AM 

1, P0O,P01 , T 00, T01 , AM ACM, AG,RG 

2 « AMASS 1 , AMASSO, T “ 

3 ,S00, CMAS, CMOM, CENR 

A -rSOr bp, -c m , -vi sc i — 

COMMON ZA, ZB, UJO, DOR IN 

6 , MASRAT, R ZERO , AF AC , CFO, "HD, F, Al , A A7ERH, T F , 


EQUATIONS CHANGE, 


ALT T 
ALIT 
A L T T 
ALIT 
ALIT 


FP, GP 


7, SDLOSS , AS D 

6 , AUG1 , U1 , UCENT 

REAL MASRAT 

-DIMENSION XI 3N G0( 3t , "F < 10T " 

1 , AN ( 10 ) , DC 10 ) ,0 2N( 10 ) ,D1D( 10 ) ,D2D(10) 

, A 1 ( 10) , TF(10> 

5 ,FP( 10) ,GPC 10) 

—UR * XC 1) 

Y 2 *= X ( 2 ) 

DOR = X C 3 ) 

Z1 * 0. 



Z 3 - 0 

74 = 0 

Z 5 = 0 

— rs--~o 

' CTR * -.05 

$l 0 -tt-Iv/C W-SO0*UR*UR*Y£*Y2) 

S 20 = Y 2* S 10 

— 0S20Y2 = S10 >5 20* S 20* ?• * SCO* UR *UR 

DS20UR = S2O*S2O*Y2*S0O*UR*2. 

DO 89 I'll 10 

CTR = CTR + *1 

E S - TF ( I ) 

GO TO (60, 61,62) ,J 

i > c AN<n - Y2 >F ( M 

D2N ( I ) = 1 • 

GO TO -"63 “ 

61 AN ( I ) * ( Y 2 + F(I))*(Y2 + F(I)) 


ALIY 

ALIT 


AL IT 
ALIT 
ALIT 
ALIT 
ALIT 
ALIT 
ALIT 
ALIT 
ALIT 

ALIT 
AL IT 
ALIT 
ALIT 
ALIT 


ALIT 

ALIT 

ALIT 

ALIT 


n 

2'i 

^ ) 
A) 
5 I 


67 

7 ■} 


3 3 
70 
13 ) 
11 > 
12) 
13: 

1 A) 
150 
163 

2 70 
23 3 
7 77 

3 ) 
310 


3 A 3 
35 ) 
3 53 
370 
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D2Nm »“2.*(Y2 + Ft IM ALIT 330 

GO TO 63 At IT 390 

~62 AMH* 1Y2 + Ft I)) *( l.~ + T*FS! - At IT 4 J ) 

D2N t I ) 3 1 • + T*FS ALIT 41 ) 

“6^3 DU) • »T* F 5-SOO* Uft Y ? +F1 IH*tY2*F< II f ALIT 4 ?> 

010( I) = -2.*SOO*UR*(Y2 + Ft I ) ) * ( Y 2 ♦ F ( I ) 1 ALIT 430 

D2Dm~- -2.*SOO*UR*UR*(Y2 + F( I) ) ALIT 44 ) 

Zl * Z1 + 0.2*AN( I I /DC I ) AL IT 

i2 = 72 -o«2*A^m*DiDm/(nm*Dmt altt 


\ *CTR 

Z3 = Z3 + 0.2*D2Nt M/Dt I) 

1 * CTR 

Z4 = 14 .0.2^AN(n*020tH/(D(II*0m) 

1— *CTft ^ 

89 CONTINUE ALIT 490 

D THREE-EOUAT-ICNS- AND-THEIi^ -DERIVATIVES WITH RESPECT- TO UR, Y2r OOR - ALIT 5>} 


RHOO = POO / ( RG* TOO I 

-c POO LOSS AT INLET OR NOT 

P02 = POO 

P02 = POO - 9# 16* ( AM AS 50* 32* 2 /^. 3T2) **2.*.072/ l 32. 2*RH00 ) 

AGG = AG/ ( AG- I* ) AL IT 51 ) 

-S3- - -T>- SOO*UR*UR*Y2*Y2 AL I t~ 5- V) 

Y4 * S3*+AGG*P02/P00 

DY4UR ---Y4*AGG/$3*2,*SOO*UR*V2*Y2 

DY4Y2 = -Y4*AGG/S3* 2.*SwO*UR*UP *Y2 

GO TO ( 2 3? 24 ♦ 24 1 ) # J — — AL I T *5 5 ) 

23 FFtJI = ( COR*COR*U°* ( 71-S2CI ♦ UR*S20)*Y4 ALIT 570 

jGGM )-^(OOR+OOR*t 71—520 ) — * — S 203 *Y 4 ALIT -530 


l +(DOR*nnR*UR*IZ?-r)S20UR) + UR *P S20UR ) * Y4 

2+(0QR*DQR*UP*{ ZI-S20 ) + UR* S 20 ) *P Y4UR 

GGt 2 ) = ( DOR*DOR*UR* ( Z3 + 74-0S20Y2 J * UR*0S20 Y2) * Y4 

2+<DOR*DOR*UR*( 71-S20 ) + UR*S 20) *f>Y4Y2 

G G I 3 I = 2* *PCR*UR*t Z 1- $20 ) *Y4 

DIVON * l* /CMAS 

GG TO 88 

-24 FF( J ) = { f>OR*DOR*UR*UR*< 7 1-S20*Y2 > 

1 ♦ UR*UR*Y2*S20)*V4 + IY4-1. ) *. 5/ < AGG+SOO > *AFAC 

G Gin - ( OO R* DOR *2 o *UR*( 7l-S?0*Y2l 

1 + D0P*D0R*UR*UR*tZ2-DS20UR*Y2) 

2- ~ 2.*UR^5 20*Y2-v-H^R*UR*DS20UR^Y-2 = 3+Y4 

3 + DY4UR*. 5/f AGG*500)*AFAC 

4 + ( DOR*DOR *UR *UR *( 21-5 20 *Y2) * UP*UR* Y2*S20 ) *D Y4UR 

GGt 2 ) =( DOR*OOR*UR*UR *( 73 + Z4-0S20 Y2*Y2 -520 ) 

— I + UR*UR*< Y2*0S 20 v 2 +520) )*Y4 — — 

3 + DY4Y2*. 5/( AGG* 500 ) *A F A C 

4 + (DOR*OOR*UR+UR*( Z1-S20*Y2) + U R*UR *Y2* S20 ) *D Y4 Y2 
GGt 3 ) - 2. *COR*UR *UR* ( 7 1 -S20*Y2)*Y4 

OIVDN = 1. /C M CM 

GG TO 88 

241 0V1 = 1. - $0O*UR*UR*Y2*Y2 

WV = SOO*TOC/TO! 

DV2 = 1. -WV*UR*UR*(1. + Y2)*f l"i"'»~Y2 V - — * 

HV2 = 1./DV2 

FF(J) = CV1*0V2 

GGt 1 ) = -HV2*SOO*2. *UR*Y2*Y2 

- ■— 1 ♦DV1*DV2*DV2*WV*2.*UR*(1.+Y2)*(1.+Y2> ~ 

GGt 2 ) = -DV2*SOO*UR*UR*2.*Y2 

1 -♦DVl*0V2*DV2*WV*UR*UR*2.*l 1.+Y21 — 

GGt 3 ) = 0 . 


ALIT 53) 

ALIT 4- JO 

ALIT 61) 

— ALIT 6 2) 

ALIT 6)) 


ALIT 64) 


ALIT 


At 

IT 

740 

AL 

tt 

770 

AL 

IT 

730 

AL 

IT 

7 TO 

AL 

TT 

300 

AL 

IT 

310 

AL 

IT 

320 

AL 

IT 

330 

AL 

IT 

340 

AL 

IT 

35) 
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_p C0NST= (poi/P02>**( ( AG- 1 • 1 /Af^l 

DIVDN = l./PCONST 

88 CONTINUE 

GO TO ( 25 , 26t 27 ) , J 

25 FFCl) ~ FFfl t -CMAS 

GO TO 28 

26 FFl 2) * FFI2I - CMCH 

GO TO 28 

27 CONTINUE 

F F( 3 ) = FF ( 3 ) -PCONST 

“28 CONTINUE 

C TO NORMALIZE THE FUNCTIONS R F { ) AND ITS DERIVATIVES GG( > 

FF ( J ) * FF ( J )*D I VDN 

DO 111 1=1,3 

“GGCTI * GGm*OIVDN 

111 CONTINUE 

C CALCUt ATI ON" OF TNT Tt AL 'VALUE OF "DELTA/ R 1 APPROX I MATE I 

DORTN * SORT U CMOM *0.5*Y2*Y2 - Y2* S20) / ( Z 1 -S2 0*Y2 >) 


ALIT 863 
ALIT 37) 
ALIT 330 
ALIT 391 
ALIT 9); 
ALIT 11 ) 
ALIT 9 20 
ALIT 940 
ALIT 950 


— RETURN " 

END 

SUBROUTINE SUB ( XX , R , DR , J , R ZERO) " 

DIMENSION 4(50), X(50), D{50), DDC50) ALITD60 

GO-TO-f 2*2tl-lTJ “ ALI T 13 to 

1 1=0 ALTT1330 

C“ — READ X AND D I AM ■ ~ ALIT 109 3 

C XX=X/RZERO ,R=R/RZEPO ALIT11K. 

WRITE (6, 5) • ” " 

5 FORMAT ( 15X ,12HPR0EILE DATA , / / , 1 OX , 5HX ( FT ) , 9X ,7HD I A { F T ) ,/) 

— 11 1 = I ♦ t — ALITlll ) 

READ (5, 10) XII), At I ) ALTTU2; 

— 10 FORMAT (2F15.4) 

IF (All)) 11,12,11 AUTUVl 

12 !ENO=I^l — -- ALITH5) 

SCALE = 1 • /RZERO 

[30 50 1=1, TEND * At IT 1 1 7) 

WRITE(6,10)Xm,Am 

All )= A( I )*« 5*SCALE - ALITin J 

50 xm = X{ d*scale alituto 

WR I T£ I 6 , 1 5 I — 

15 FORMAT!///) 

— I^T- — — — — At IT l 2 ) > 

C ALT T 121 ) 

C FIND R AND DR ALIT12 20 

C ALIT1231 

2 IMXX-X(U) 20,23, 22 ' 4LIT1240 

20 1=1-1 AL IT 1 2 50 

tf-1 XX-X ( I) )”20 f 23, 23 A At I T 1 269 

23 I B= I AL IT 1 271 

1A- 1 ♦ l At IT 1239 

GO TO 24 ALIT12D 

-22 I=TM At IT 1 3 39 

IF (XX-X(l)l 25,23,22 ALIT1319 

25 14=1 ALIT132) 

IB=I-1 At IT 1 3 39 

— 24 DR= ( A( I A)-A( IB) ) /( X| I A )-XU8)) — ALIT l 341 

R=AI IB)+(XX-X( IB))* DR ALIT1351 

RETURN — ALTT136) 

END A L I T 1 3 7 J 

StJBROUT TNE Rt JNGE~ (Nt~X 7 “Yy~TOt , “YMI N,“H, XOUT, MARK) ALIT 390 

C FIRST ORDER DI FF* EO. ROUT INE — A0JUS1S STEP SIZE ALIT 11 0 
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230 

250 

10 — 

20 

30 


DIMENSION YI 10) ,YMIN110 ) , TOL( 10 > » SUR( 10) , XHUTI 100) ,MARK( 5) 

-- YA,50, » FAI50I, FB ( 50 ) * FC(50»* YKEEP<50> 

KBIG = i 

— 

NCOUNT = 1 5 

J = MARK ( 1 ) ... 

MAX = M ARK ( 2 ( 

DO 250 I -*-l, N 

SUB(I) = TOL (1 1/32.0 

IF IN AX-= J )— 20, 30 ,—30 


RETURN 

A XOUT ( J ) - X 

B = ABS (2. F-6*XJ 
IF t A ♦ B) AO, 35, 35 
35 IF (A - B) 50, 50, 60 
40 J * J t 


-50 


GO TO 10 
CONTINUE 
CALL PRINT 
J *“J ♦ 1 


(N, XOUT, Y, OY, J) 


GO TO 10 

-60 1 F C A 1 .5 *H> — 70,— 70 ,80- 

70 H « A 

GO TO 1000 


IF IA - 3. *H) 
H-» v5*A 


80 
-40“ 

C 

C 

lOC-fc — XA -» X ♦ H/3. 

XB = X ♦ . 5*H 


90, 1000, 1000 


DO RU N G E-KtfTY E-MG R SON I NT5CR ATI PN 


CALL CIFFEO (N, X, 
X * X ♦ H 

-00 -1030 — * — t* — N — 


OV ) 


YKEFP(I) = Y ( I ) 

FA ( I ) * H*DY( H 

1030 YAII) = Y( I ) + FAlD/3. 

CALL OIFFEO (N, XA, YA, OY1 

DO 1040 I « 1, N 

1040 YA( t)=Y ( I ) -*- FA ( I t /6, H*DYf! I /6. 

CALL DIFFEO (N, XA, YA , OY) 

DC 1050 — I-t -l, N 

FBIII = H* DY ( I ) 

1C5C -Y A I I ) Y { I ) -K . 1 25*FA ( I I— ♦ — r 575*FBM ) 
CALL DIFFEO (N,XR, YA, DY) 

~0G 1060 1 — » “1 N — - 


-1060 


FC( I > = H*DY( I ) 

YA(IV Y( n + «5*FA{I) !• 5*FB f TO— #*- -2-« *FC f f) 

CALL DIFFEO (N.X, YA, DY) 

00 1130 1 = 1, N 


Y(I) = Y(I) + FAtn/6, * • 666666667+FC ( II + H*DY(I)/6. 
-1C 61 — U * YU) 


I F ( ABS C U 1 - YMINC I I) 

1090 KLOW = 2 

F * • 2* AB S (U - YAC I> ) 

IF ( E AB$ ( T QL ( I )*U) 

1100 KBIG * 2 

GC TO 1130 

1110 IF ( F - ABS C SUB C I )*U) 


1130, 1130, 1090 


> — 1 M 0 .- 1110 • 1100 


) 1130, 1120, 1120 


ALT T 

ALTT 33) 
ALIT 340 
ALIT 35) 
—ALIT "3 SO 
ALIT 37) 
ALIT 33) 
ALTT 39J 
ALIT 40) 
ALIT 410 
“ At IT 4 20 
ALIT 430 

ALIT 440 

ALIT 450 

ALIT 460 

ALIT 470 

A tr l T 4 3 0 

ALIT 490 

ALIT 500 

ALIT 51) 
- ALIT 520 
ALIT 530 

At IT — 54 0 

ALIT 550 

ALTT 560 

ALIT 570 

ALIT 530 

ALIT 593 

At IT- 60^0 

ALIT 313 

ALTT 620 

ALIT 63) 
- ALIT 640 
ALIT 65) 
ALIT- 660 

ALIT 670 

ALIT 630 

ALIT 690 

ALIT 700 

ALIT 710 
—At IT - 7 20 
ALIT 733 

“ALIT 740 

ALIT 750 

ALIT 760 

ALIT 770 

ALIT 79) 

ALIT 3)0 

ALIT 31) 

ALIT 920 

ALIT 930 

ALIT- 340 

ALIT 950 

ALIT 96) 

ALIT 3 7 j 

ALIT 390 

ALIT 99) 
ALIT 900 
ALIT 910 
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1120 KBTWN =2 

1130 CONTINUE _ 

GO TO ( 100, 1135), KLOW 

1135 GC TO (1180, 1140), KBIG 

1140 NCOUNT = NCOUNT -~t 

IF (NCOUNT) 1150, 1150, 1170^ 

•1150 PRINT 1160. X, H ' " 

PRINT 1165, (I, YU), DY(I), I = 1, N) 

1160 FORMAT (58H4STEP SIZE HALVED 15 TIMES CONSECUTIVELY SINCE LAST 

jujj / PROGRAM TERMINATED AT X ■= E16.8, BH, H * , El 6 

2//3H l , 13X, 4HY( I ) , 1 6 X , 5HDYII),//) 

1165 FORMAT ( I 3 , 7X , E 16. 8 , 4X , F 1 6. 8 ) ' ’ 

1170 KBIG = 1 

IF (H - B) 1176, 1172, 1172 

1172 X = X - H 

H * »5*H 

DO 1174 1=1, N 

1174 Y( I) = YKEEPl I ) 

KBTWN = 1 __ 

KLOW = 1 

GO TO 1000 

1176 M = 15 - NCCUNT - 
PRINT 117e, M, X, 


H 


1176 


118C 

1 IPO 
120C 

C 

C — 

c 


PRINT 1165, (I, Y ( I ) , DY ( I ) , I = l, N) 

RETURN 

FORMAT (41H4STEP SIZE BECAME TOO SMALL FOR COMP UTER . /2FH IT HAS 
1EN HALVED , 12, 21H TI”ES CONSECUTIVELY. /29H PROGRAM TERMINATED 

-2 x j— £16. B, 8H, - H = , E16.-fif/V3H I, 13X, 4HYI.I), 16X, 

35HDY ( I) , // ) 

NCOUNT = 15 ~ 

GO TO (11P0, 1200), KBTWN 

H = 2.*H — - 

KBTWN =1 

~ KLOW-^-1 


CHECK FOR INTERMEDIATE PRINT OUT 
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140 


ICO XOUT ( 50 ) =X 

JK=0 

C C ALL“PR I NT ( N, XOUT , Y ,DY , JK ) 

IF (JK) 10,10,20 . . , 

FORMAT (5H X = , E16.B, 4X, 4HH- , E16.B, 1 IX , 1HI , 13X, 4HY( I 

116X, 5HDY ( I)/) 

FORMAT (55X, 13, 7X, 2IE16.8, 4X)) 

END _ 

SUBROUTINE -SIM0(A,B,N,KS1 

DIMENSION A ( 1 ) , B ( 1 ) 

C ...»..«« .«»« »»»»»» ••••••••••• 

c 

c 

t 


SUBROUTINE SIMO 


C 

e 

c 

c- 

c 

c 


PURPOSE 

- OBTAIN SOLUTION OF A SET T)F SIMULTANEOUS LINEAR EQUATIONS, 
AX = B 


USAGE 

CAtt — S IMQ( A tO »NfX^SI 


ALIT 923 
ALIT 933 
ALIT 943 
ALIT 95) 
ALIT 959 
ALIT 97) 
ALIT 933 
ALIT 99) 
ALITU)0 
PPIALITI )1 ) 
9, ALIT132) 
ALIU )3 J 
ALIT1 >49 
ALIT1353 
AL!Tl)6) 
ALIT13 7 ) 

■ — ALI T 1930 
ALIT1 ) 9 3 
At IT 1 1)3 
ALT Till 3 
ALIU12) 
AL I T1 1 3 ) 
ALIT 114) 
ALITH53 
ALIT116) 
ALIT1173 
REALIT11VJ 
ATALIT1 1 93 
ALIT12 )‘J 
ALIT121 ) 
ALIT122 ) 
AL I T 1 233 
ALIT 124) 
AL l Ti 25 3 
ALTT1250 
ALIT1273 
ALIT1230 
ALIT1293 
AL I T 1 339 
ALIT131 ) 
A L IT 1 3 ? O 
AL I T 1 3 1 3 
\ t AL IT 1 349 
ALIT135) 
ALIT 136=3 
AL IT 1 370 
-■ ALIU350 
ALIT137 ) 
AL I T 1 3^) 
►.•ALIT 1393 
ALIT1433 
ALIT1410 
ALIU42) 
AL I T 1430 
ALIT 1449 
AL IT 1450 
AL IT 1 46) 
AL IT 147 ) 
ALTT1439 
AL I T 1490 
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oJjouuJjuOuuu6oiiuouuoi> uouoououuo 


C OESCR IPT ION OF PARAMETERS - - - - 

C A - MATRIX OF COEFFICIENTS STORED COLUMNWISE. THFSE ARF 

t DESTROYEO-IN THE COMPUTATION. THE SIZE OF MATRIX A IS 

C N BY N, 

C B- VECTOR OF ORIGINAL CONSTANTS (LENGTH N>. THESE ARE 

c REPLACED BY FINAL SOLUTION VALUES, VECTOR X. 

C ' tf - _ NUMBER OF EOUA T IONS AND VARIABLES. N MOST RE • GT. ONF. 

C KS - OUTPUT DIGIT 

C £ FOR A NORMAL SOLUTION 

I FOR A SINGULAR SET OF EQUATIONS 


ALIT150 ) 
ALIT151 ) 
ALITI5>3 
ALIT153J 
ALTT15V3 
ALIT1550 
ALITISSJ 
ALIT1570 
AL IT 1 5 3 3 
ALITiSRO 


REMARKS 

MATRIX A MUST BE GENERAL* ■ 

IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS^ 


AN ALTFRNAT IVE SOLUTION MAY 3F OBT A I NED BY USING MATRIX 
INVERSION (MINV) AND MATRIX PRODUCT (GMPRO). 


SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE— 


METHOD- 


METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PIVOTAL 
“01 VT SOP ir EACH ST AGE~ OF EL TM I NAT ION CONS! S TS OF INTEPCHANG I 
ROWS WHEN NECESSARY TO AVOID DIVISION BY ZERO CR SMALL 

“ELEMENTS. 

THE FORWARD SOLUTION TO OBTAIN VARIABLE N IS DONE IN 
N "ST AGES • THE BACK SOLUTION FOR THE OTHER VARIABLES IS 
CALCULATED BY SUCCESSIVE SUBSTITUTIONS, FINAL SOLUTION 
~VALUES*ARE ^DEVFtrOPEfT IN VECTORS, WITH VARIABLE 1 IN B(l) , 
VARIABLE 2 IN B(2), ........ v VARIABLE N IN B(N). 

IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF 0.0, 

THE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. THIS 
TOLERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT 


-- ALITlsn 
AL IT IS 10 
ALIT1620 
ALITI63 J 
ALITISA ) 
ALIT165 ) 

ALTTtSS ) 

AL IT1S7 G 
ALTTISBJ 
ALITlSVj 

- ALITI7 JO 

alititlo 
NGAL T t 1 7 > ) 
AL TT 17 30 
ALTT1740 
ALTT17S ) 
AL IT17SU 
ALIT1770 

- ALIT173C 
AL I f l 7 00 
ALIT13)3 
ALTT131J 

. AL I TIB ?0 
AL I T 1 3 3 ) 

> 6 AL I T l 3 St) 


FORWARD SOLUTION 


T0L=0 .0 

KS -0 

JJ=-N 

..... oo 6 5 j-* ^ M 

JY = J + 1 

— JJ=JJ*N+1 .... 

B1GA=0 

TT- J J -J — — — — — — 

DO 30 I=J,N 

C SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN 



I J= I T + 1 

IF( AflS( BIGA)-ABS( A( IJ ) ) I 20, 30, 30 

20 BIGA = A{ IJ ) 

IMAX=I • 

30 CONTINUE 

C TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX) 

IFUBS(BIGA’)-TOL) 35, 35, SO 


Atm 350 
AL IT 1330 
A L I T 1 3 >0 
ALITIDJu 
AL I T 1 0 1 ) 

AL I r 1° 30 

\L t T 1 9 3C 
ALIT19S0 
AL IT195 ) 
A L 1 T 1 9 Si) 
AL IT 197 > 

AL I T 1^3) 

ALIT199") 
ALIT7) ) 0 
ALIT2U ) 
- AL I T 70 2u 

A L I T 7 } 3 0 

ALIT70SU 

AL IT 20 5 ) 
- ALIT2)6, 
ALIT237 ) 
AL IT 7 ) 30 
AL I T 7 ) 90 

AL I T 210 ./ 

AL T T 21 1 ) 
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35 KS-1 

RETURN 

INTERCHANGE ROWS IF NECESSARY 

40 I1*J*N*(J-2J 

f T* IM AX- J 

00 50 K = J f N 

T 1* 1 1+N ‘ 

12-1 1*TT 
SAVE=A( 1 1 ) 

A ( 1 1 ) = A f T 2 ) _ „ • 

A( 12 ) =S AVE 

0IV1DE EQUATION BY LEADING COEEFICIENT 

— 5C“Ai ITT^AT I X V/BITjA ' " “ 

S AVE= 0( I MA X ) 

B(INAXt=B( J) 

R ( J ) * SAVE/ 3 I C A 

ELIMINATE NEXT VARIABL E __ 

IMJ-N) 55,70,55 

55 10S*N*(J-I1 

OC 65 I X = J Y , N 

I XJ= I QS + 1 X " 

I T = J- I X 

00 6^“ JX- J Y, N 

IXJX*N*( JX-ll+IX 

JjX=IXJX+IT 

60 A( IXJX)xA( IXJX)-( A( IXJ)*A(JJX)I 

— 65 BUX)=&(IX)-( 8 {J)*AUXJU 

C 

C— ^ACK SOLUTION 

C 

— 70 NY=N-1 

I T = N*N 

OC 80 J * 1 , NY “ * 

l A* I T- J 

TB=N-J 

I C-N 

— DO 80 K=l,J 

B ( I B ) = B ( IB ) -A ( I A ) +8 ( I C J 

I A= ! A-N 

8C IC=IC-1 

RETURN * — - 

SUBROUTINE VB04A (M,N,F,X,E,A,E$CALE, I PR I NT , MA XKf G , T A , W ) 

DIMENSION A ( 3t 10 ) , W(10), F(3), E(3)t X(3) 

— — ic=o - 

I KFG=0 

I P= I PRINT 

IP5ET=-IPRINT*l IPRINT-l ) 

NN=N+N 

NP=N+ 1 

NNP = N^NP 

I CON= l 

1 FS - 1 

I S=2 


ALIT21 20 

ALU213 j 

ALIT214 ) 

ALIT2150 

ALIT2160 

A L I T 2 1 7 ) 
AL IT 21 BO 
ALIT21? ) 
AL IT 23 30 
, AL IT 22 1 i 
ALIT222 ) 

A L I T 2 2 3 0 

AL IT 2 24 ) 

AL IT 2250 
ALIT226 i 
ALIT22TO 

AL T T 2 28 3 

AL I T 2270 
ALIT23D) 
ALIT231 3 
ALIT232J 
ALU233 j 

ALIT23V3 

ALIT2353 
AL IT 7 3 5 3 
ALTT2373 
ALIT2333 
ALIT230 ) 
AL I T 24 ) J 
AL I T 24 1 3 
ALIT242) 
A L t T 24 3 ) 
AL IT 2440 
M IT ?45 ) 
AL I T 24 5 3 
A L I T 2 4 7 3 
AL I T 2 + 3 ) 
ALIT2473 
ALIT253 ) 
AL I T 2 5 1 3 
AL I T 2 5 20 
AL I T 2 53 ) 
ALIT2543 
ALIT255 3 
AL IT 2 5 6 ) 
ALI T 757 3 
AL I T 2530 

AL IT 20 3 3 
ALIT21 ) ) 
AL I T 21 1 0 
ALIT2123 

ALIT213) 

A L I T 7 1 4 ) 
ALTT2153 
ALIT216 ) 
ALIT21 70 
ALIT2130 

AL IT?!' 7 ) 

ALIT22 JO 
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— I KK=N 

DO 3 1=1, N 

KK = KK + 1 

DO 2 J=1,I 

—2 - AM •JV=0. 

3 W(KK>=0. 

F F=0 • 

I KFG= IKFG+1 

DO 4 K=l,M 

CALL CALXFG (K,X,F,W> 

KK=N 

DO 5 1 = 1, N 

KK=KK+1 

DO 6 J*l # I 

— A(I, J1 = A( I.J)+W(T >*W(J> 

5 W(KK)=W(KK)-W(t)*Fm 

4 FF=FF+F(K)*F<K ) 

GO TO ( 7,8) , IS 

— B DC 45 I =2 » N 

IK = I - 1 

DC 45 J = I,N 

45 AUK, J) = M J, TK) 

CAtL' SFNISTTA,ir,M, TK) 

KK=NP 

DO 9 1 = 1, N — 

W(I)=W(KK) 

9 KK=KK*1 — 

EM*0. 

OO-lO-T^ir** — 

W (KK ) =0 • 

DO 11 J = l » N "■ 

11 W (KK ) =W ( KK )+A( I »J ) * W ( J ) 

^ r MAX1( EM,ABS(W(KK)/Et I))) 

10 KK=KK+ 1 

1 F-( EM- E SC ALE 1-1^, 1 2, 13 


ALIT7710 
AL TT 222 ) 
ALIT2233 
AL I T 2 2 40 
AL IT >?50 
ALTT226) 

A L I T 2 2 70 

ALIT2230 
AL I T 2200 
AL IT2330 

A L IT 2 310 

ALIT2330 

ALU7TT0 

ALIT2340 

ALIT2350 

ALTT2353 
ALIT2370 
ALIT2333 
AL IT 2 300 
' ALIT240J 
ALIT2410 
AL IT742 ) 

ALIT2430 

ALIT244) 
AL I T 2 4 5 0 
AL IT 2450 
ALIT2470 
AL I T 24 30 

At IT 249 0 

ALIT2510 
ALIT2S1 ) 
AL IT 2 5 2 0 
ALU253 ) 
A L I T 2 540 

AL IT2 5 50 


13 EM=ESCALE/EM 

KK =NNP “ ~~ 

DO 14 1=1, N 

W (KK ) =W ( KK )*EM 

14 KK=KK+ 1 

WTO -15— ~~ 

12 IF (EM-1.) 16,16,15 

16 GO TO (27,18), ICON • 

27 I CON = 2 

GO TO 1 7 ' 

18 IF (IPRINT) 19,19,20 

— 19 RETURN 

20 WRITE (6,21 ) 

—21 FORMAT (//5X,18HV804A FINAL VALUES) 

I FS= 2 

GO TO 22 

15 l CON= 1 

— 17 1 Fi I KF 0-M AXK FG ) 28 ,29 ,29 — 

29 WR IT E ( 6 , 30) IK F G 

30 FORMAT ( / / 5 X , 5HVB04 A » I 5 , 1 7HM CALLS OF CALXFG) 

GO TO 18 

—28 IP = IP*!PRINT 

IF (IP) 23,23,22 

— 22 WRITE! 6,26 )IC, IKf G 

26 FORMAT ( //5X, 9HITFR ATION, 14, 1 10 , 1 7HM CALLS OF CALXFG) 


ALU256 ) 
AL1T2570 
A L I T 2 5 3 ) 
AL IT 2 59 } 
ALIT2400 
A LIT 2610 
ALIT2420 
AL I T 26 3 ) 
AL I T 3 6 40 
A L I T 2 6 50 
A L I T 2 6 6 0 
- AL IT 2570 
ALIT2630 
AL IT 26 70 
AL I T ? 7 0 0 
AL IT 27 1 3 
AL IT 27 20 
AL IT 27 30 
AL IT 2 740 
ALIT2750 
AL T r 2 7 6 ) 
ALIT2770 
ALTT273 j 
AL IT 2790 
ALIT2300 
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- ^wm(6,24 ) ( X { I) , ! - 1 f N > 

24 FORMAT ( /5X,9HVARTA3LFS/, (5E24* 14)) 

- WRITE(6,25HF( I), I=l,M> 

25 FORMAT { /5X , 9HFUNCT IONS / , ( 5E24. 14 ) > 

GO TO-(^8» 191 , IG$ - 

38 I P= I PS ET 

23 TC-I C+i 

ITEST=3 

— ffx = FF — 

XP=0. 

X ocu 

1S=1 

- 7 GG*0 . 

KK=NNP 

00 31 I *NP , KN 

GG=GG-W( I)*W(KK) 

-21— KK-KK+1 

GC=GG*GG 

GALL VD02A C ITEST,XC ,FF , GG, 6,0. ,0. 3, 1,4- 

GO TO (32, 33,33,33, 33), ITEST 
-32 XP*XC-XP 


v 4Lir2siO 

ALIT2333 

ALIT2333 

ALIT2340 

ALIT2359 

ALIT2360 

A L IT 28 70 

ALT T 2330 

ALIT2390 

4LIT29J > 

-At! T 291 > 

ALIT292) 

*LIT2939 

AL IT 2940 

AL IT 29 5t> 

AL I T 2950 

At I T 29 70 

ALIT2930 

nALIT2990 

ALU3D33 
AL I T 30 10 


KK=NNP 

00 34 1*1, N 

X(I )=X( I ) *XP*W ( KK ) 

-34 KK=KK+1 

XP = XC 

* GC TO 1— J 

33 IF (FF-FFX) 39,40,40 

-39-I FtABS( EM^XC) -^1, )36r3t>,~3T 

27 I CON= 1 

~26 IF (XP-XC) 35,8,35 

25 I S= 2 

GO TO 3 2 

40 WRTTE(6,41 ) 

-41 FORMAT ( //SX,33FVft04A- ACCURACY “CANNOT E* AC M IF V E 01 

GC TO 18 

— END ; — — 

SUBROUTINE VD02A (I TEST , X ,F , G , MAX FUN, ABS ACC , REL ACC , XSTEP ) 

ABSF ( XY Z)=AB$(XY/I 

S IGNF( ABC, XY7 ) = S I GN ( ABC , X Y7 > 

— — SQRTfM XYZ)=SQRT(XYZ f 

GO TO ( 1,2,2) , ITEST 

— 2 I S=I TEST 

ITEST *1 

1 INOl 

XINC=XSTEP 

— — MOl 

GC TO (4,4,10) , IS 

— 3 MC^HOl ” 

IF (MC-MAXFUN) 4,4,5 

5 I TEST = 4 ■ 

11 F=FA 

— X*X A — - 

G=GA 

— 4 RETURN 

10 I S=2 

MC=0 — 

1 IF (G) 6,7,6 

T - ITEST =5" 

GO TO 4 

—6 GO TO (8, 9), IS 

9 XA=£ 


A LIT 30 29 
A L IT 30 3 ) 
A L I T 3940 
AL IT 30 5 0 
AL IT 3 0 60 
ALTT3070 
ALTT303) 
ALT T 3090 
ALTT31 10 
ALIT3110 
At I T 31 2 ) 
ALITU30 
ALIT3149 
- ALIT315) 
ALIT314) 

ALTT3170 

ALT T 3130 
AL I T 3 1 9.) 
AL I T 320 ) 

AL 1732 19 

ALTT3220 
AL IT 32 30 
AL IT 3 240 
ALTT3250 
ALIT3263 

ALIT3270 

A L I T 3 2 30 
AL IT320 ) 
ALTT330 0 
ALIT3310 
AL IT3320 

AL I T 3 3 3 ) 

AL I T 3 340 
AL I T 3350 
ALIT3360 
AL IT 3379 
AL 17 3 3 39 
AL IT3390 
AL IT340 ) 
ALTT3419 
AL I T 34 2 0 
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GA *G 

1 S= 1 

12 X = X A— S I GNF { X INC ,GA ) 

— X 1 NC= X I NC ♦ X I NC 

GO TO 3 

8 IF i F-F A) 13,13,14 

13 OUM=FA 

FA-F ” 

F = DUM 

— OUM^GA 

G A = G 

— G = OUM 

OUM=XA 

XA*X 

X=OUM 

14 IF (GAM X-XA1 1 15,16,16 

15 I INC = 2 

X INC=X 

16 Z = 3#*(FA-F ) / ( X- XA ) +G+GA 

W^Z^Z-G*GA ; 

IF <W) 20,20,17 

17 W-S IGNF (SQRTFfW) - , X-XA) 

XP=X-(X-XA)*(G+W-Z) M G-GA+W+W) 

IF ((XP-XA)*GA) 18, 19,20 " 

18 GO TO < 21, 221 , ITNC 

51 IF ( ABSFCXP-XAl-ABSMXTNC )) 23, 23 , 12 
22 IF ( AB5F( XP-XA)-APS F( XINC-XP ) ) 23,23,24 

24 X=0.5*(X!NC*XAI ~ 

IF ( (X-XA) *{XINC-X) ) 25,25,3 

25 TTEST-3 

GC TO 11 

"23 X^XP 

IF ( ABSF( XP-X A I-ABSF ( ABSACC) ) 19, 19,26 

26 IF ( ABSFt XP-XA ) ~ARS F t XP*R EL ACC ) 1 19,19,3 

19 ITE$T=2 

GO TO 11 

20 GC TO (12,241,1 INC 

TO 

SUB90UT TNE S PN I ST ( U , I , J , X ) 

DIMENSION U(3, 31 ,Vt 3,3) 

I F ( I - 3) 2,1,2 

— 2 I = J ♦ K ^ 

WR I T E ( 6 , 4 ) 

— 4 T r CRMATT14HO"MOOIFT"-V8C4A— 1 — 

STOP 

— 1 Vtifl t — = — 0 C 2, 2)*U C 3 , 3 > ”** (Jf 3y2!**2 
V ( 2 , 2 ) * UCl, 11*U(3,31 - U<3, 11**2 
V( 3 , 3 ) U(1,1)*UI2,2) - U< 2, 11**2 

V ( 2 , 1 1 * U(3,1)*U(3,2) ~ U(2,U*UI3,3) 

V ( 1,-2 ) — = V ( 2 , 1 1 

V ( 3 , 1 ) = 11(2,11*11(3,21 - U(2,21*U (3,1 I 

vc 1,-31-*- v ( 3 , i j ““ 

V ( 3, 2 ) = U(3,l)*U(2t 11 - U( 1 t i I *U ( 3,2 1 

V(2,31 * V(3,21 — “ 

OET * U(l,l)*V(l,l) + lf(2,l)*V(2, 1) «■ U ( 3 ,1 1 *V ( 3 , 1 I 

0C-5-1r-»“1 ,3 

DO 5 M= 1 , 3 

5 U(l,M) » V(L,M1/DET “ 

RETURN 

END 


A L IT 34 30 
AL IT344C 

A L IT 34 50 

AL I T 3460 
ALIT3470 
AL I T 3 4 30 
AL IT 3490 
ALIT 3500 
■ - AL T T 3510 
AL I T 3 5 20 
AL I r 3 5 3 0 
AL I F 3 54 0 
ALIT3550 
AL I T 3 560 

A L T T 3 5 70 

AL IT 35 30 
AL I T 3 590 
AL I T 3 6 00 
ALIT3410 
A L I T 3620 
ALIT3630 
A L IT 3 540 
AL IT 36 50 
AL l T 3560 
ALIM670 
AL IT 3 6 30 
ALIT3690 
ALIT3700 
ALI T 3710 
A L I T 3 7 20 
ALIT 3730 
AL 1^3740 
AL I y 3 7 50 
ALIT376G 
\Lt r 3^70 
AL I T 3 7 3 0 
ALIT379Q 
ALIT 3300 
ALIT 38 10 
ALIT3320 
A L I T 3 8 3 0 
AL I T3840 
ALI73850 
AL I r 3 B 60 
ALIT3370 
A L I T 3g 80 
A L I T 3 8 90 
AL I r 3 900 
ALIT3910 
ALT r 392Q 

AL ! T 39 30 

ALIT 3940 
AL l 7 3 950 
AL I T 3 960 
ALIT3970 
AL I T 3980 
ALI T 3990 
ALIT4000 
ALIT 40 10 
ALIT4020 
AL I T 40 30 
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B. 4 Typical Sets of Input and Output Data 
Input Data: 


10 

.9950,9850. 9 7 50. 9650.9550,9450. 9350.9 25 0,91 50. 9C c 0 


.33 

.21 



.000 22 3 

1.400000 .000160 .223000 



M = 

21.0 



. 500 

5.2981 147.0 1 1 

2. 50 

, 00046 

14.7 552. 

0 348.0 1267.0 .113 2.375 

1.40 

53. 2 

.0 

.4^5 



.6315 

.3035 



.3815 

.3037 



. 76 86 

.4462 



. 2148 

.4463 



. 2981 

. 5796 



.000000 

0.000000 




10 

. 99 50.9 850.97 50. 9650.9550. 94 5 0. 9 350 . 9250. 91 50. 9050 


.33 

.000223 

M= 

.21 

1.400000 .000160 

17.0 -.02 

. 223000 



. 500 

5. 2981 147.0 

2 

1 

2.00 

. 00046 

1 4.7 552. 

0 348.0 1267.0 . 

113 

1.885 

1.40 

53. 2 
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TABLE 2 

PRESSURE TAP LOCATIONS AND FINAL MIXING TUBE DIMENSIONS 


Static 
Pressure 
Tap No. 

Stagnation 

Pressure 

Traverse 

No. 

Pressure Tap 
Location 
figure 4 
x-inches 

Dimensionless 

Location 

x/Ro 

(R 0 = 2. 670") 

1 


0.46 

0.172 

2 


2.21 

0.828 

3 


4.71 

1.76 

4 

1 

6.71 

2.51 

5 


9. 71 

3.63 

6 

2 

12.21 

4.57 

7 


14.71 

5.51 

8 

3 

17.21 

6.45 

9 


20.46 

7.66 

10 


22.21 

8.32 

11 

4, 5 

24.71 

9.25 

12 


27.21 

10.19 

13 


29.71 

11.13 

14 


32.21 

12.06 

15 


34.71 

13.0 

16 


37.21 

13.94 

17 


39.71 

14.87 

18 


42.21 

15.81 

19 


44.46 

16.65 

20 

6 

47.21 

17.68 

21 


49. 71 

18.62 

22 


51.69 

19.36 

23 


57.69 

21.60 

24 


63.69 

23.85 

25 


69.69 

26.09 


Measured Mixing 
Tube Dimensions 


HU 

— 

wTEm 

x/R 

o 

0 

5.341 


19.578 

3.643 

7.34 

28.578 

3.645 


45.224 

5.355 

16.9 

50.578 

5.356 

18.95 

63.578 

6. 956 

23.8 
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TABLE 3 

STATIC PRESSURE VALUES MEASURED ALONG THE MIXING TUBE 
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VELOCITY AND TEMPERATURE PROFILES AT TRAVERSE STATION 1 
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Figure 1 Jet Mixing in a Converging-Diverging Duct 



Figure 2. Velocity Profiles in the Transition Zone 
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Figure 4 Mixing Tube Geometry 



Figure 5 Jet Pump Test Rig 




Primary Flow Measurement 
t irifice 
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Figure 7 Mixing Tube Static Pressure Variations 




Local Radius Ratio 


Figure 8 
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Measured Velocity Profiles in Mixing Tube 
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Comparison of Analytical and Experimental Mixing Tube Static 
Pressure Variations ( x c nrf , * 2. 5, entrainment reduced by 2%) 
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-Centerline Temperature - R U -Centerline Velocity - ft/sec 
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Figure 14 Comparison of Analytical and Experimental Values of 

Centerline Temperature and Velocity Along the Mixing Tube 
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Velocity - fps 


Test Data 



Y~ - Local Radius Ratio 

Figure 15 Comparison of Experimental and Analytical 
Velocity Profiles 
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